Hello !
I’m currently an intern trying to code in dolfinx 0.8.0 (python 3.11.13) to create a library of solvers, and tools to facilitate usage on complex meshes. I’ve been having a hard time finding working datasheet or API documentation that are working in FenicsX as tutorials are not always updated.
As I am an intermediate user of FenicsX, I’m tackling a full-coupled thermoelastic solver using mixed element spaces. There are probably a lot to optimize, especially regarding the syntax.
I’m getting heavily inspired by the Thermo-elastic evolution problem (full coupling) — Computational Mechanics Numerical Tours with FEniCSx tutorial, made by Jeremy Bleyer.
I solved all the error messages and the solver runs, however I am now stuck with a divergent solution where spotting the error is a real torment. I’ve verified the definitions of boundary condition in order to narrow the potential errors, but I can’t find any in the solver part itself.
I’m then linking my code in order that you could help me spot (and potentially solve) the problems.
Here is my simple mesh and a function to locate elements (to assign BC) :
#Dimensions en m (X,Y,Z) et éléments du maillage (X_M,Y_M,Z_M)
X,X_M=5,5
Y,Y_M=1,1
Z,Z_M=1,1
domain = mesh.create_box(MPI.COMM_WORLD, [np.array([0, 0, 0]), np.array([X, Y, Z])],
[X_M, Y_M, Z_M], cell_type=mesh.CellType.hexahedron)
#repérage des facettes (surfaces à l'air libre) se fait directement dans locate_entities_boundary, donc on lui donne un array rempli de True.
def boundary(x) :
return np.full(x.shape[1], fill_value=True,dtype=bool)
dim=domain.topology.dim-1 #dimension de surface
#Stockage de toutes les facettes extérieures
all = mesh.locate_entities(domain, dim, boundary)
markers = []
values = []
#Fonctions de détections de surfaces, via les coordonnées (X,Y,Z) avec tolérance
def left(x): return np.isclose(x[0], 0,atol=1e-8)
def right(x): return np.isclose(x[0], X,atol=1e-8)
def front(x): return np.isclose(x[1], 0,atol=1e-8)
def back(x): return np.isclose(x[1], Y,atol=1e-8)
def bottom(x): return np.isclose(x[2], 0,atol=1e-8)
def top(x): return np.isclose(x[2], Z,atol=1e-8)
#Stockage de chaque ENTITE dans facets avec correspondance boundary_conditions dans values_all.
boundary_conditions = [(left, 1), (right, 2), (front, 3), (back, 4), (bottom, 5), (top, 6)]
for func, marker in boundary_conditions:
facets = mesh.locate_entities_boundary(domain, dim, func)
markers.append(facets)
values.append(np.full(len(facets), marker, np.int32))
facets_all = np.hstack(markers)
values_all = np.hstack(values)
mt = mesh.meshtags(domain, dim, facets_all, values_all)
#ATTENTION on localise les entités (=éléments) et non les noeuds crééons donc un lien
domain.topology.create_connectivity(dim, domain.topology.dim)
domain.topology.create_connectivity(domain.topology.dim, dim)
#La fonction dirichletbc prends les dofs en entrée, d'où la fonction suivante
def entities_to_dofs(Vsub, domain, dim, entities):
dofs =fem.locate_dofs_topological((Vsub, Vsub.collapse()[0]), dim, entities)
return dofs
Here are the function spaces and associated functions :
# Definition des espaces adapté à chaque grandeur physique
Vu = basix.ufl.element( "P", domain.basix_cell(), 1, shape=(3,) ) #le shape donne la dimension, "P"2 donne le type d'élément et basix_cell repère les éléments.
Vt = basix.ufl.element("P", domain.basix_cell(), 1) #Idem avec la température
#Définition de l'espace fonctionnel mixte
V = fem.functionspace(domain, basix.ufl.mixed_element([Vu, Vt]))
#Cofonction servant à l'extraction des fonctions mixtes
v = fem.Function(V)
#Sous-fonctions de chaque sous-espaces fonctionnels
u , T = v.sub(0), v.sub(1)
u.name , T.name = "Déplacement" , "Température"
#COPIE indépendante des sous-espaces pour manipulation de données
U , _ = V.sub(0).collapse()
V_t , _ = V.sub(1).collapse()
u_func = fem.Function(U)
T_func = fem.Function(V_t)
# Température de référence (CL sur T)
Tref_value=293.15
Tref=fem.Function(V_t)
Tref.x.array[:]=Tref_value
dt = fem.Constant(domain, 1.0) #incrément de temps
Here are the BC definition and variables :
#Stockage des DL pour visualisation
dofs_T=[]
dofs_u=[]
#Fixed
uu = fem.Function(U)
uu.x.array[:] = 0
#dofs
CLAMPED=entities_to_dofs(V.sub(0), domain, dim, mt.find(1))
CLAMPED_ID=CLAMPED[1]
dofs_u.append(CLAMPED_ID)
bc_CLAMP=fem.dirichletbc(uu , CLAMPED ,U)
#Traction
u_lim=fem.Function(U)
for i in range(U.dofmap.index_map.size_local):
u_lim.x.array[3*i] = 0.1 # x seulement
#dofs
right_dofs=entities_to_dofs(V.sub(0), domain, dim, mt.find(2)) #Tableau des indices des éléments adjacents à la surface droite
right_dofs_ID=right_dofs[1]
dofs_u.append(right_dofs_ID)
bc_U=fem.dirichletbc(u_lim, right_dofs, V.sub(0))
#Ext. temp.
uT=fem.Function(V_t)
uT.x.array[:] = Tref_value + 10
#dofs
combined_dofs = dolfinx.fem.locate_dofs_topological((V.sub(1), V_t), dim, all)
bc_T=fem.dirichletbc(uT, combined_dofs , V_t)
#Uniquement pour le stockage
combined_dofs_local = dolfinx.fem.locate_dofs_topological( V_t, dim, all)
dofs_T.append(combined_dofs_local)
# used for post-processing
bottom_T_dofs = fem.locate_dofs_topological(
(V.sub(1),V_t), dim, mt.find(5))
bc=[bc_CLAMP,
bc_T,
bc_U ]
E = fem.Constant(domain, 70e9) #Mpa
nu = fem.Constant(domain, 0.3) #sans unité
rho = fem.Constant(domain, 2700.0) #kg/m3
k = fem.Constant(domain, 0.)
alpha = fem.Constant(domain, 0.)
cV = fem.Constant(domain, 910.0) #J/K/kg
gdim=domain.topology.dim
lmbda = E * nu / (1 + nu) / (1 - 2 * nu)
mu = E / 2 / (1 + nu)
kappa = alpha * (3 * lmbda + 2 * mu)
I made an totally independant part to visualize the boundary conditions with pyvista. It is not mandatory :
# --- CL en déplacement ---
u_vis = fem.Function(V.sub(0).collapse()[0])
u_vis.x.array[:] = 0.0 # par défaut : pas de CL
#Bloc de vérification des dimensions (on évite de comparer les DL de température dans le lot en utilisant les sous-espaces)
n_dofs = u_vis.x.array.size
all_dofs = np.unique(np.hstack(dofs_u)) #concatène tous les dofs où des CL de déplacement s'appliquent
assert all_dofs.max() < n_dofs, "indices non-locaux"
#Définition du champ local
bs = 3
node_ids = np.unique(all_dofs // bs) #comme les dofs spatiaux se suivent par 3, la division euclidienne les trie de manière unique
arr = u_vis.x.array.reshape(-1, bs) #créé artificiellement une correspondance 3-uplet (dofs) scalaire (noeuds)
arr[:] = 0.0 #initialise les valeurs nodales à 0
arr[node_ids, :] = 1.0 #attribution de la valeur 1 à chaque noeuds précontraint
u_nodal = arr.mean(axis=1) #Passage à forme nodale (avec moyenne si CL sur un dof et pas 3)
# --- CL en température ---
T_vis = fem.Function(V.sub(1).collapse()[0])
T_vis.x.array[:] = 0.0
node_idT=dofs_T #Dofs des CL sur T
ar=T_vis.x.array
ar[:] = 0.0 #initialise les valeurs nodales à 0
ar[node_idT] = 1.0 #attribution de la valeur 1 à chaque noeuds précontraint
T_nodal = ar #valeurs nodales=valeurs dofs ssi dim=1
#Conversion vers PyVista pour affichage
# Déplacements : surface 0/1
topo_u, cell_u, geom_u = plot.vtk_mesh(u_vis.function_space)
grid_u = pyvista.UnstructuredGrid(topo_u, cell_u, geom_u)
grid_u["BC_displacement"] = u_nodal
# Température : surface 0/1
topo_T, cell_T, geom_T = plot.vtk_mesh(T_vis.function_space)
grid_T = pyvista.UnstructuredGrid(topo_T, cell_T, geom_T)
grid_T["BC_temperature"] = T_vis.x.array.real
# Affichage côte à côte (PyVista)
p = pyvista.Plotter(shape=(1, 2))
p.subplot(0, 0)
p.add_text("BC Déplacement", font_size=10)
p.add_mesh(grid_u, scalars="BC_displacement", cmap="plasma", show_edges=True)
p.add_axes()
p.subplot(0, 1)
p.add_text("BC Température", font_size=10)
p.add_mesh(grid_T, scalars="BC_temperature", cmap="coolwarm", show_edges=True)
p.add_axes()
p.link_views()
p.show()
And here is probably where lies the error, the variational formulation :
v = fem.Function(V) #fonction de formulation faible
u, Theta = ufl.split(v) #en (déplacement,température) d'où le split.
v_ = ufl.TestFunction(V) #fonction test pour la résolution
u_, Theta_ = ufl.split(v_)
dv = ufl.TrialFunction(V) #fonction d'incrément
V_aux = fem.functionspace(domain, ("DG", 1))
s_old = fem.Function(V_aux, name="Previous_entropy") #Servira dans l'itération temporelle (suite arithmétique)
v.sub(1).interpolate(uT)
v.sub(0).x.array[:] = 0.0
def eps(u):
return ufl.sym(ufl.grad(u))
sig = (
lmbda * ufl.tr(eps(u)) * ufl.Identity(gdim)
+ mu * eps(u)
- kappa * Theta * ufl.Identity(gdim)
)
j = -k * ufl.grad(Theta)
s = cV / Tref * Theta + kappa / rho * ufl.tr(eps(u))
s_expr = fem.Expression(s, V_aux.element.interpolation_points())
s_old.interpolate(s_expr)
mech_res = ufl.inner(sig, eps(u_)) * ufl.dx
therm_res = (
rho * Tref_value * (s - s_old) / dt * Theta_ - ufl.dot(j, ufl.grad(Theta_))
) * ufl.dx
Res = mech_res + therm_res
J = ufl.derivative(Res, v, dv)
#forme matricielle de la Jacobienne
Jac=fem.form(J)
A = fem.petsc.create_matrix(Jac)
#Forme vectorielle de l'erreur résiduelle
Resid=fem.form(Res)
L = fem.petsc.create_vector(Resid)
The definition of my solver :
solver = PETSc.KSP().create(domain.comm)
solver.setOperators(A)
solver.setType("gmres")
solver.getPC().setType("ilu")
du = fem.Function(V)
And finally the Newton iteration loop :
i = 0
max_iterations = 25
v.sub(0).interpolate(u_func)
v.sub(1).interpolate(T_func)
bs = U.dofmap.index_map_bs # blocksize (devrait être 3)
#initialisation de v :
P=(X_M+1)*(Y_M+1)*(Z_M+1)
v.x.array[:] = 0
v.x.array[-P:] = Tref_value
print("INITIALISATION\t","v=" ,v.x.array)
while i < max_iterations:
# Assemble Jacobienne et erreur résiduelle
with L.localForm() as loc_L:
loc_L.set(0)
A.zeroEntries()
dolfinx.fem.petsc.assemble_matrix(A, Jac, bcs=bc)
A.assemble()
dolfinx.fem.petsc.assemble_vector(L, Resid)
dolfinx.fem.petsc.apply_lifting(L, [Jac], [bc], x0=[v.x.petsc_vec])
L.ghostUpdate(addv=PETSc.InsertMode.ADD_VALUES, mode=PETSc.ScatterMode.REVERSE)
#Boundary condition :
dolfinx.fem.petsc.set_bc(L, bc, v.x.petsc_vec)
# Scale l'erreur par -1
L.scale(-1)
L.ghostUpdate(addv=PETSc.InsertMode.INSERT_VALUES, mode=PETSc.ScatterMode.FORWARD)
# Résolution du problème
solver.solve(L, du.x.petsc_vec)
du.x.scatter_forward()
##Visualisation
# 1) formes & tailles
v_test = v.x.array
print("shapes: v", v_test.shape)
# 2) premières valeurs (contraste)
print("v[-4:]:", np.array(v_test)[-4:])
print("v[4:]",np.array(v_test)[:4])
# Itération
v.x.array[:] += du.x.array
i += 1
u_vals = u_func.x.array.reshape((-1, bs)) # shape (ndofs_u, 3)
T_vals = T_func.x.array
# Compute norm of update
correction_norm = du.x.petsc_vec.norm(0)
print(f"\tIteration {i}: Correction norm {correction_norm}")
norm_T_2 = np.linalg.norm(v_test)
norm_T_inf = np.max(np.abs(v_test))
print("Normes mixtes (2,inf) :", norm_T_2, norm_T_inf)
if correction_norm < 1e-10:
break
I’ve removed a lot of parts where I added prints and visualizations of different variables. These were only for me to find the cause of divergence. My code goes to the max iterations really fast with a huge divergence on the correction norm, giving negative values into my solution vector too.
v[-4:]: [-4.28897665e+48 4.54580087e+46 -4.54690223e+46 4.28896650e+48]
v[4:] [-2.47568404e+46 8.87744845e+46 8.69098051e+46 -2.46270201e+46]
Iteration 25: Correction norm 1.6076882919887018e+51
Mixed norms (2,inf) : 5.801950325714887e+50 3.5551309979642035e+50
I wonder if the cause doesn’t come from the mixed element spaces themself and their specific syntax.
Sorry for the gargantuan code I provided, as I said it still lacks brevity.
Hopefully, someone can lead me to the real problem.
Thank you in advance !