Hello all,
I am develoing a code to solve a multiphysics poblem and blow is the set of equation I am solving in each element.
I am using a mixed elment formulation and also using lagrange multipliare to force a constraint in my problem.i.e. 
.
my code works if I use a direct solver for my problem, but I need to use very fine mesh, so I need to use an iterative solver. Based on my search, using block solver is a good way to use in this kind of problem, like the on in fenics22-tutorial.
So I made changes in my code but I got the following error:
Here is the code I am running:
from mpi4py import MPI
from petsc4py import PETSc
import numpy as np
import pyvista
import math
from dolfinx.fem import (Constant, Function, functionspace, assemble_scalar, dirichletbc, form, locate_dofs_geometrical,
locate_dofs_topological, assemble, locate_dofs_topological, extract_function_spaces, bcs_by_block)
from dolfinx.fem.petsc import assemble_matrix, assemble_vector, apply_lifting, create_vector, set_bc, LinearProblem, assemble_matrix_nest, assemble_vector_nest, apply_lifting_nest, set_bc_nest
from dolfinx.io import XDMFFile, VTXWriter, gmshio
from dolfinx.mesh import create_box, locate_entities, meshtags, locate_entities_boundary, CellType
from petsc4py.PETSc import ScalarType
import meshio
from ufl import (FacetNormal, Identity, TestFunctions, TrialFunctions,
div, dot, dx, inner, lhs, nabla_grad, rhs, sym, indices, as_tensor, Measure, split, action,
SpatialCoordinate)
from basix.ufl import element, mixed_element
from dolfinx import default_scalar_type
import time
import dolfinx
print(
f"DOLFINx version: {dolfinx.__version__} based on GIT commit: {dolfinx.git_commit_hash} of https://github.com/FEniCS/dolfinx/")
start = time.time()
Fl = 1e-8
Fs = 1e-5
Fp = 1e-4
# ....................
C_star = Fl ** 2 / Fs
h_star = 1 / Fl
mu_star = Fp / Fs
perm_star = Fp ** 2 / Fs
q_star = Fp / (Fs * Fl)
mo_star = Fp / (Fl ** 2)
diff_star = 1 / (Fl ** 2)
p_star = Fl ** 3
F_star = Fs
qe = (1.6021 * 1e-19) * q_star
# -------------------------Material Properties
E = 1 * 152e9 * C_star
nu = 0.33
la = E * nu / ((1 + nu) * (1 - 2 * nu))
mu = E / (2 * (1 + nu))
perm = 146.6e-9 * perm_star
epsilon = 11e-9
non_local = 10e-9 / Fl
mu11 = 0.5 * 121e-6 * mu_star
mu12 = 0.5 * 121e-6 * mu_star
mu44 = 0.0 * mu_star
# ......................................... Reading mesh file
H = 0.76e-3 / Fl
mesh= create_box(MPI.COMM_WORLD, [np.array([0, 0, 0]), np.array([H,H, H])],
[20, 20, 20], cell_type=CellType.hexahedron)
# ......................................... Define Loading
Force = 200.00000 / Fs
Area2 = 1.2769000000000007e-06 / Fl ** 2
tr_S2 = Constant(mesh, ScalarType((0, 0, -Force / Area2)))
noder = FacetNormal(mesh)
# NN = Constant(mesh, ScalarType((1.0, 1.0, 1.0)))
# ......................................... Define Elements
v = element("Lagrange", mesh.basix_cell(), 1, shape=(mesh.geometry.dim,)) # Displacement Element
s = element("Lagrange", mesh.basix_cell(), 1, shape=(3, 3)) # Relaxed strain Element
lm = element("Lagrange", mesh.basix_cell(), 1, shape=(3, 3)) # Lagrange multipliers Element
p = element("Lagrange", mesh.basix_cell(), 1) # Electric Potential Element
# ......................................... Define mixed element
MIXEL = mixed_element([v, p, s, lm])
Space = functionspace(mesh, MIXEL)
# ......................................... Define Trail and Test Functions
(u, p, gradu, lamu) = TrialFunctions(Space)
(del_u, del_p, del_gradu, del_lamu) = TestFunctions(Space)
#########################################################3
# ....................................... Mark boundary regions
tol = 1e-16
def Bottom(x):
return np.isclose(x[2], 0, atol=tol)
def Top(x):
return np.isclose(x[2], H, atol=tol)
# ....................................... ds
boundaries = [(1, lambda x: np.isclose(x[2], 0.0, atol=tol)), (2, lambda x: np.isclose(x[2], H, atol=tol))]
facet_indices, facet_markers = [], []
fdim = mesh.topology.dim - 1
for (marker, locator) in boundaries:
facets = locate_entities(mesh, 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 = meshtags(mesh, fdim, facet_indices[sorted_facets], facet_markers[sorted_facets])
ds = Measure("ds", domain=mesh, subdomain_data=facet_tag)
# ....................................... BCs
S0, S0_to_Space = Space.sub(0).collapse()
# Clamp BC for left side
clamp = Function(S0)
clamp.x.array[:] = 0
clamp_facets = locate_entities_boundary(mesh, mesh.topology.dim - 1, Bottom)
# .........Fixed BC
clamp_dofs = locate_dofs_topological((Space.sub(0), S0), mesh.topology.dim - 1, clamp_facets)
bc0 = dirichletbc(clamp, clamp_dofs, Space.sub(0))
# # .........Ground BC
P0, P0_to_Space = Space.sub(1).collapse()
ground = Function(P0)
ground.x.array[:] = 0
ground_facets = locate_entities_boundary(mesh, mesh.topology.dim - 1, Bottom)
ground_dofs = locate_dofs_topological((Space.sub(1), P0), mesh.topology.dim - 1, ground_facets)
bc1 = dirichletbc(ground, ground_dofs, Space.sub(1))
# ........................
bc = [bc0, bc1]
# ......................................... Define tensor
delta = Identity(len(u))
i, j, k, l, m, n = indices(6)
# Strain Tensor
def StrainT(u):
return as_tensor((1. / 2. * (u[i].dx(j) + u[j].dx(i))), [i, j])
# Stress Tensor
def Elas(la, mu):
return as_tensor(la * delta[i, j] * delta[k, l] + 2. * mu * delta[i, k] * delta[j, l],
(i, j, k, l)) # Elasticity tensor
def StressT(u):
return as_tensor((Elas(la, mu)[i, j, k, l] * StrainT(u)[k, l]), [i, j])
def RStrianGrad(gradu):
return as_tensor((1. / 2. * (gradu[k, j].dx(i) + gradu[j, k].dx(i))), [i, j, k])
def SGElas(la, non_local, mu):
return as_tensor((la * (non_local ** 2) * delta[i, j] * delta[l, m] * delta[k, n] +
2 * mu * (non_local ** 2) * delta[i, l] * delta[j, m] * delta[k, n]),
(i, j, k, l, m, n)) # Strain Gradient Elasticity Tensor
def HStress(gradu):
return as_tensor(
(SGElas(la, non_local, mu)[i, j, k, l, m, n] * RStrianGrad(gradu)[l, m, n]),
(i, j, k))
# Flexoelectricity Tensor
fl1 = np.zeros((3, 3, 3, 3))
fl1[0, 0, 0, 0] = 1
fl1[1, 1, 1, 1] = 1
fl1[2, 2, 2, 2] = 1
flexo1 = Constant(mesh, PETSc.ScalarType(fl1))
flexo = as_tensor(((mu11 - mu12 - 2 * mu44) * flexo1[i, j, k, l] + mu12 * delta[i, j] * delta[
k, l] + mu12 * delta[i, k] * delta[j, l] + mu12 * delta[i, l] * delta[j, k]), (i, j, k, l))
# PREPARE FORMS
kuu = StressT(u)[j, k] * StrainT(del_u)[j, k] * dx
kua = -lamu[j, k] * del_u[k].dx(j) * dx
kpp = -perm * delta[i, j] * (-1 * nabla_grad(p))[j] * (-1 * nabla_grad(del_p)[i]) * dx
kps = -(flexo[i, j, k, l] * RStrianGrad(gradu)[j, k, l]) * (-1 * nabla_grad(del_p)[i]) * dx
ksp = -1 * flexo[l, i, j, k] * (-1 * nabla_grad(p))[l] * RStrianGrad(del_gradu)[i, j, k] * dx
kss = HStress(gradu)[i, j, k] * RStrianGrad(del_gradu)[i, j, k] * dx
ksa = lamu[j, k] * del_gradu[j, k] * dx
kau = del_lamu[j, k] * u[k].dx(j) * dx
kas = del_lamu[j, k] * gradu[j, k] * dx
Zer = Constant(mesh, ScalarType(0.0))
kaa = Zer * lamu[j, k] * del_lamu[j, k] * dx
Zer2 = Constant(mesh, ScalarType((0, 0, 0)))
fu = tr_S2[i] * del_u[i] * ds(2)
fp = inner(Constant(mesh, 0.), del_p) * dx
loadsa = as_tensor([[0, 0, 0], [0, 0, 0], [0, 0, 0]])
fs = Zer2[i] * del_gradu[i, j] * noder[j] * ds
fa = Zer2[i] * del_lamu[i, j] * noder[j] * ds
# .........................................CREATE BILINEAR AND LINEAR FORMS
a = form([[kuu, None, None, kua], [None, kpp, kps, None], [None, ksp, kss, ksa], [kau, None, kas, kaa]])
L = form([fu, fp, fs, fa])
# .........................................PREPARE PRECONDITIONER
paa = lamu[j, k] * del_lamu[j, k] * dx
PR = form([[kuu, None, None, None], [None, kpp, None, None], [None, None, kss, None], [None, None, None, paa]])
# APPLY BCs
A = assemble_matrix_nest(a, bcs=bc)
A.assemble()
b = assemble_vector_nest(L)
apply_lifting_nest(b, a, bcs=bc)
for b_sub in b.getNestSubVecs():
b_sub.ghostUpdate(addv=PETSc.InsertMode.ADD,
mode=PETSc.ScatterMode.REVERSE)
bcspace = extract_function_spaces(L)
bcs0 = bcs_by_block(bcspace, bc)
set_bc_nest(b, bcs0)
P = assemble_matrix_nest(PR, bc)
P.assemble()
# ......................................CREATE SOLVER
ksp = PETSc.KSP().create(mesh.comm)
ksp.setOperators(A, P)
ksp.setType("minres")
ksp.setTolerances(rtol=1e-9)
ksp.getPC().setType("fieldsplit")
ksp.getPC().setFieldSplitType(PETSc.PC.CompositeType.ADDITIVE)
nested_IS = P.getNestISs()
ksp.getPC().setFieldSplitIS(("u", nested_IS[0][0]),
("p", nested_IS[0][1]),
("gradu", nested_IS[0][2]),
("lamu", nested_IS[0][3]))
# Set the preconditioners for each block
ksp_u, ksp_p, ksp_gradu, ksp_lamu = ksp.getPC().getFieldSplitSubKSP()
ksp_u.setType("preonly")
ksp_u.getPC().setType("gamg")
ksp_p.setType("preonly")
ksp_p.getPC().setType("jacobi")
ksp_gradu.setType("preonly")
ksp_gradu.getPC().setType("gamg")
ksp_lamu.setType("preonly")
ksp_lamu.getPC().setType("gamg")
# Monitor the convergence of the KSP
ksp.setFromOptions()
### SOLVING
Result = Function(Space)
ksp.solve(b, Result)
assert ksp.getConvergedReason() > 0
# ......................................... Export Results
Displ, potent, DisGrad, LaMA = Result.split()
Displ.x.scatter_forward()
potent.x.scatter_forward()
DisGrad.x.scatter_forward()
LaMA.x.scatter_forward()
#...........................
uf = Displ.collapse()
pot = potent.collapse()
with VTXWriter(mesh.comm, ("Dis.bp"), [uf]) as vtx:
vtx.write(0.0)
with VTXWriter(mesh.comm, ("Pot.bp"), [pot]) as vtx:
vtx.write(0.0)
Any idea what might be wrong here?