# FEniCSx result doesn't coincide with previously known solution

Hi all. I’m new to FEniCS and FEniCSx. After many days of trying and a lot of questions here, I could write a full working code for a test problem and find a result. Unfortunately, I already know the expected solution, and the result I got is not correct. With these values of the parameters, the solution is a uniform vector field pointing to the right. However, the result looks like this:

The yellow vectors are the boundary conditions. It’s difficult to see there, but there are small vectors in the rightmost columns (I guess those more to the left have a modulus ~0) whose direction inverts from column to column. Moreover, the result doesn’t change at all with the values of the physical parameters I use or the tolerance of the solver, which converges suspiciously fast.

There must be something fundamentally wrong… I hope you or anyone can help me with this. Thanks in advance!

The code is:

``````import numpy as np
import ufl
from dolfinx import mesh, fem, nls, log, plot
from mpi4py import MPI
from petsc4py.PETSc import ScalarType

dim = 3 # dimension of the OP
a = -1; c = 0.01 # thermotropic parameters
d = 1 # stiffness

Rbeta = 1  # radius of the base (sphere)

domain = mesh.create_rectangle(MPI.COMM_WORLD, np.array([[0,0],[Rbeta,np.pi]]),[10,int(10*np.pi)])
x = ufl.SpatialCoordinate(domain)
tdim = domain.topology.dim; fdim = tdim - 1
element = ufl.VectorElement("CG", domain.ufl_cell(), 1, dim)
V = fem.FunctionSpace(domain, element) # space of vector functions with dim components

m0 = [1,0,0] # preferred value of the OP in the boundary
m00 = [1,0,0] # preferred value of the OP in the center

def surface(x):
"""Identifies the elements in the surface of the sphere."""
return np.isclose(x,Rbeta)

def origin(x):
"""Identifies the elements in the center of the sphere, a.k.a. the origin."""
return np.isclose(x,0)

BCs = []
boundary_facets = mesh.locate_entities_boundary(domain, fdim, surface) # facets in the surface
for i in range(dim):
boundary_dofs_x = fem.locate_dofs_topological(V.sub(i), fdim, boundary_facets) # i-th comp. of the surface dofs
BCs.append(fem.dirichletbc(ScalarType(m0[i]), boundary_dofs_x, V.sub(i))) # preferred value in the boundary
boundary_facets = mesh.locate_entities_boundary(domain, fdim, origin) # facets in the origin
for i in range(dim):
boundary_dofs_x = fem.locate_dofs_topological(V.sub(i), fdim, boundary_facets) # i-th comp. of the origin dofs
BCs.append(fem.dirichletbc(ScalarType(m00[i]), boundary_dofs_x, V.sub(i))) # preferred value in the origin

def Jacob(x):
"""Jacobian of the transformation from Cartesian to spherical coordinates."""
xr, xq = x
return xr**2*ufl.sin(xq)

def curl_sph(A, x):
"""Curl of a vector field depending only on r and theta in spherical coordinates."""
xr, xq = x; jr, jq = 0, 1
Ar, Aq, Af = A
curl_r = (ufl.Dx(Af*ufl.sin(xq),jq))/(xr*ufl.sin(xq)) #(ufl.Dx(Af*ufl.sin(xq),jq)-ufl.Dx(Aq,jf))/(xr*ufl.sin(xq))
curl_q = (-ufl.Dx(xr*Af,jr))/xr #(ufl.Dx(Ar,jf)/ufl.sin(xq)-ufl.Dx(xr*Af,jr))/xr
curl_f = (ufl.Dx(xr*Aq,jr)-ufl.Dx(Ar,jq))/xr
return ufl.as_vector([curl_r, curl_q, curl_f])

def div_sph(A, x):
"""Divergence of a vector field depending only on r and theta in spherical coordinates."""
xr, xq = x; jr, jq = 0, 1
Ar, Aq, Af = A
div_r = ufl.Dx(xr**2*Ar,jr)/xr**2
div_q = ufl.Dx(Aq*ufl.sin(xq),jq)/(xr*ufl.sin(xq))
div_f = 0 #ufl.Dx(Af,jf)/(xr*ufl.sin(xq))
return div_r + div_q + div_f

def fB(u, x):
"""Bulk energy density considering the thermotropic contribution and elastic distortions."""
fT = a/2.*ufl.dot(u,u) + c/4.*ufl.dot(u,u)**2
fD = d*( (div_sph(u,x))**2 + (ufl.dot(u,curl_sph(u,x)))**2 + (ufl.cross(u,curl_sph(u,x)))**2 )
return fT+fD

m = fem.Function(V) # dim-dimensional vector OP - what we are trying to find
dm = ufl.TrialFunction(V) # variation of the OP
phi = ufl.TestFunction(V) # test function

F = ufl.dot(m,m)*Jacob(x)*ufl.dx # total bulk energy fB(m,x)
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, bcs=BCs, J=JF)

solver = nls.petsc.NewtonSolver(MPI.COMM_WORLD, problem)
solver.convergence_criterion = "incremental"
solver.rtol = 1e-6
solver.report = True

log.set_log_level(log.LogLevel.INFO)
n, converged = solver.solve(m)
assert(converged)
print(f"Number of interations: {n:d}")
``````

The output is simply:

``````Number of interations: 2
2022-08-04 16:04:57.760 (   0.124s) [main            ]              petsc.cpp:677   INFO| PETSc Krylov solver starting to solve system.
2022-08-04 16:04:57.767 (   0.131s) [main            ]              petsc.cpp:677   INFO| PETSc Krylov solver starting to solve system.
2022-08-04 16:04:57.770 (   0.134s) [main            ]       NewtonSolver.cpp:36    INFO| Newton iteration 2: r (abs) = 5.10179e-16 (tol = 1e-10) r (rel) = 6.20363e-17(tol = 1e-06)
2022-08-04 16:04:57.770 (   0.134s) [main            ]       NewtonSolver.cpp:255   INFO| Newton solver finished in 2 iterations and 2 linear solver iterations.
``````

Does anyone have any suggestion on how to start looking for the problem? Does any part of the code look weird? I find it very suspicious that the solver converges so fast, but I couldn’t find any information on this on the internet.

I would suggest starting by trying to use the “incremental” method for the newton solver, and set the underlying ksp solver to use a direct solver.

Hi @dokken, thanks for your reply. The incremental method is already in use. As for the ksp solver, could you please explain how to pass this option? I added the line `solver.ksp_type = "preonly"` (apparently this is the way to call a direct solver?), but nothing different happened. I also tried other ways of passing the option (`solver.set('ksp_type', 'preonly')` and
`solver.PETScOptions.set('ksp_type', 'preonly')`) but I got diverse error messages.

Sorry, i meant using the “residual” method.

With the residual method and that `solver.ksp_type = "preonly"` line it converges even faster, and the plot looks exactly like the one in the first post.

``````Number of interations: 1
2022-08-09 09:54:07.746 (   7.239s) [main            ]       NewtonSolver.cpp:36    INFO| Newton iteration 0: r (abs) = 8.00001 (tol = 1e-10) r (rel) = inf(tol = 1e-06)
2022-08-09 09:54:07.748 (   7.242s) [main            ]              petsc.cpp:677   INFO| PETSc Krylov solver starting to solve system.
2022-08-09 09:54:07.756 (   7.249s) [main            ]       NewtonSolver.cpp:36    INFO| Newton iteration 1: r (abs) = 2.48777e-18 (tol = 1e-10) r (rel) = 3.02507e-19(tol = 1e-06)
2022-08-09 09:54:07.756 (   7.250s) [main            ]       NewtonSolver.cpp:255   INFO| Newton solver finished in 1 iterations and 1 linear solver iterations.
``````

This would indicate that the problem you have created above is actually linear,
which one can quite easily see from the variational form:

This means that `dF=2*ufl.dot(m, phi)*Jacob(x)*ufl.dx` and that you can use a Linear solver to get the solution of dF.

1 Like

Thanks @dokken! The simplified energy I had used as an example lead indeed to a linear problem. When I put the real expression, everything worked fine. I guess simpler is not always better 