Hi everyone,
I am currently comparing the performance of direct and iterative solvers in my code, but I am encountering an issue where the iterative solver does not produce the same results as the direct solver.
Has anyone experienced a similar issue or have any suggestions on what might be causing the discrepancy?
Here is my code with iterative solver:
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)
from dolfinx.fem.petsc import assemble_matrix, assemble_vector, apply_lifting, create_vector, set_bc, LinearProblem
from dolfinx.io import XDMFFile, VTXWriter, gmshio
from dolfinx.mesh import create_unit_square, locate_entities, meshtags, locate_entities_boundary, create_box, 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
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
# ......................................... Parameters
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
H = 0.76e-3 / Fl
mesh = create_box(MPI.COMM_WORLD, [np.array([0, 0, 0]), np.array([h, h, h])],
[10, 10, 10], cell_type=CellType.hexahedron)
# ......................................... Define Loading
Force = 200.00000 / Fs
Area1 = 7.398400000000025e-06 / Fl ** 2 # Bottom
Area2 = 1.2769000000000007e-06 / Fl ** 2
tr_S1 = Constant(mesh, ScalarType((0, 0, Force / Area1))) # Bottom
tr_S2 = Constant(mesh, ScalarType((0, 0, -Force / Area2)))
n = FacetNormal(mesh)
NN = Constant(mesh, ScalarType((1.0, 1.0, 1.0)))
# ......................................... Define Elements
v = element("Lagrange", mesh.basix_cell(), 2, shape=(mesh.geometry.dim,)) # Displacement Element
s = element("Lagrange", mesh.basix_cell(), 1, shape=(3, 3)) #
lm = element("Lagrange", mesh.basix_cell(), 1, shape=(3, 3))
p = element("Lagrange", mesh.basix_cell(), 1) #
# ......................................... Define mixed element
MIXEL = mixed_element([v, s, lm, p])
Space = functionspace(mesh, MIXEL)
# ......................................... Define Trail and Test Functions
(u, gradu, lamu, p) = TrialFunctions(Space)
(del_u, del_gradu, del_lamu, del_p) = TestFunctions(Space)
# ....................................... Mark boundary regions
tol = 1e-16
def Bottom(x):
return np.isclose(x[2], 0, 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(3).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(3), P0), mesh.topology.dim - 1, ground_facets)
bc1 = dirichletbc(ground, ground_dofs, Space.sub(3))
# ........................
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])
VarForm = ((la * delta[i, j] * delta[k, l] + 2. * mu * delta[i, k] * delta[j, l]) * StrainT(u)[k, l]) * StrainT(del_u)[i, j] * dx - lamu[j, k] * del_u[k].dx(j) * dx
VarForm += -perm * delta[i, j] * (-1 * nabla_grad(p))[j] * (-1 * nabla_grad(del_p)[i]) * dx
VarForm += (mu11 * gradu[j, j].dx(i) + 2 * mu12 * gradu[j, i].dx(j)) * (-1 * nabla_grad(del_p)[i]) * dx
VarForm += (-mu11 * (-1 * nabla_grad(p))[i] * delta[j, k] - 2 * mu12 * (-1 * nabla_grad(p))[j] * delta[i, k]) * (1. / 2. * (del_gradu[k, j].dx(i) + del_gradu[j, k].dx(i))) * dx
VarForm += (non_local ** 2) * ( la * gradu[l, l].dx(i) * delta[j,k] + 2* mu * gradu[j, k].dx(i)) * (1. / 2. * (del_gradu[k, j].dx(i) + del_gradu[j, k].dx(i))) * dx
VarForm += lamu[j, k] * del_gradu[j, k] * dx + del_lamu[j, k] * u[k].dx(j) * dx + del_lamu[j, k] * gradu[j, k] * dx
# ......................................... Preparing the form
a = VarForm
L = tr_S2[i] * del_u[i] * ds(2)
#.................................................SOLVER
bilinear_form = form(a)
linear_form = form(L)
A = assemble_matrix(bilinear_form, bcs=bc)
A.assemble()
b = assemble_vector(linear_form)
apply_lifting(b, [bilinear_form], [bc])
b.ghostUpdate(addv=PETSc.InsertMode.ADD, mode=PETSc.ScatterMode.REVERSE) # if run in parrallel
set_bc(b, bc)
b.ghostUpdate(addv=PETSc.InsertMode.INSERT, mode=PETSc.ScatterMode.FORWARD)
#-----------------
solver1 = PETSc.KSP().create(mesh.comm)
solver1.setOperators(A)
solver1.setType(PETSc.KSP.Type.CG)
pc1 = solver1.getPC()
pc1.setType(PETSc.PC.Type.HYPRE)
pc1.setHYPREType("boomeramg")
# --------------------
Result = Function(Space)
solver1.solve(b, Result.vector)
print(solver1.getConvergedReason())
#-----------
Displ, DisGrad, LaMA, potent = Result.split()
uf = Displ.collapse()
pot = potent.collapse()
# ......................................... Export Results
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)
end = time.time()
print('Time', end - start)
and for solving with Direct solver, I use:
# problem = LinearProblem(a, L, bcs=bc,
# petsc_options={"ksp_type": "preonly", "pc_type": "lu",
# "pc_factor_mat_solver_type": "mumps"}) #, "mat_mumps_icntl_24": 1, "mat_mumps_icntl_25": 0})
# Result = problem.solve()