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[0],Rbeta)
def origin(x):
"""Identifies the elements in the center of the sphere, a.k.a. the origin."""
return np.isclose(x[0],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.
