Hi,
I am using complex dolfinx 0.7.3 installed through coda in wsl/ubuntu. I have noticed that MUMPS solver is never able to solve a problem. Using other direct solvers such as the default LU or UMFPACK is able to solve the problem without any issues.
I am unsure if it is a problem with petsc, the conda package or the solver needs special configuration to work which I don’t know? I will appreciate if more experienced users could share their thoughts.
Below is an MWE that solves scalar Helmholtz equation in a box. The expected solution is a propagating wave through the geometry. At least on my machine switching into MUMPS never leads to a result.
import dolfinx, petsc4py, mpi4py, gmsh, ufl
from dolfinx.io import gmshio
import numpy as np
def create_geom(fov_x=0.2, fov_y=0.2, fov_z=1):
gmsh.initialize()
gmsh.option.setNumber("General.Terminal", 0) # disable output messages
gmsh.clear()
gdim = 3
model_rank = 0
fov = gmsh.model.occ.add_box(-fov_x/2, -fov_y/2, -fov_z/2, fov_x, fov_y, fov_z)
gmsh.model.occ.synchronize()
all_doms = gmsh.model.getEntities(gdim)
for j, dom in enumerate(all_doms):
gmsh.model.addPhysicalGroup(dom[0], [dom[1]], j + 1) # create the main group/node
# number all boundaries
all_edges = gmsh.model.getEntities(gdim - 1)
for j, edge in enumerate(all_edges):
gmsh.model.addPhysicalGroup(edge[0], [edge[1]], j + 1) # create the main group/node
gmsh.model.mesh.generate(gdim)
gmsh.model.mesh.refine()
model_rank = 0
mesh, ct, ft = gmshio.model_to_mesh(gmsh.model, mpi4py.MPI.COMM_WORLD, model_rank, gdim)
# gmsh.fltk.run()
gmsh.finalize()
return mesh, ct, ft
#%% input paramters
port_bnd_idx = (6, 5) # incident, transmitted
k0 = 2*np.pi/0.5
fe_degree = 3
mesh, ct, ft = create_geom()
Omega = dolfinx.fem.FunctionSpace(mesh, ("CG", fe_degree))
u = ufl.TrialFunction(Omega)
ut = ufl.TestFunction(Omega)
# bilinear form terms
F = (-ufl.inner(ufl.grad(u), ufl.grad(ut)) + k0**2*ufl.inner(u, ut))*ufl.dx
F += -ufl.inner(dolfinx.fem.Constant(mesh, petsc4py.PETSc.ScalarType(0.0)), ut)*ufl.dx # source term
#%% Port boundary at the top (excitation) and bottom surfaces
ds = ufl.Measure("ds", domain=mesh, subdomain_data=ft) # boundary differential element
n_3d = ufl.FacetNormal(mesh)
ki_3d = ufl.as_vector((0, 0, -k0)) # incident wave
kr_3d = ufl.as_vector((0, 0, k0)) # reflected wave, if any
kt_3d = ufl.as_vector((0, 0, -k0)) # transmitted wave
# incident port
bc_inc = ufl.inner(1j*kr_3d[2]*u, ut)*ds(port_bnd_idx[0])
# transmission port
bc_tr = ufl.inner(ufl.dot(n_3d, 1j*kt_3d*u), ut)*ds(port_bnd_idx[1])
F+= bc_inc + bc_tr
# excitation at the incident port
u_inc = 1
bc_exc = ufl.inner(ufl.dot(n_3d, 1j*(ki_3d-kr_3d)*u_inc), ut)*ds(port_bnd_idx[0])
F += bc_exc
a = ufl.lhs(F)
L = ufl.rhs(F)
rtol = 1e-3
petsc_options = [
{"ksp_type": "preonly", "pc_type": "lu", "monitor_convergence": True},
{"ksp_type": "preonly", "pc_type": "lu", "pc_factor_mat_solver_type":"umfpack", "monitor_convergence": True},
{"ksp_type": "preonly", "pc_type": "lu", "pc_factor_mat_solver_type":"mumps",
# "mat_mumps_icntl_2": "1",
# "mat_mumps_icntl_4": "4",
# "mat_mumps_icntl_22":"1", # out of core factorization
# "mat_mumps_icntl_24":"1", # detect null pivot rows
# "mat_mumps_icntl_28":1, # detect null pivot rows
# 'mat_mumps_cntl_3':0.01,
# "mat_mumps_icntl_1": 1,
# "mat_mumps_icntl_10":15,
# 'mat_mumps_icntl_24':1
},
]
solver_idx = 0
from dolfinx.fem.petsc import LinearProblem
Lin_system = LinearProblem(a, L, bcs = [], petsc_options=petsc_options[solver_idx])
solver = Lin_system.solver
import time
t1 = time.perf_counter()
u_sol = Lin_system.solve()
t2 = time.perf_counter()
print("Solved in {}s".format(t2-t1))
it = solver.getIterationNumber()
print("Constrained solver iterations {0:d}".format(it))
solver.view()
with dolfinx.io.VTXWriter(mpi4py.MPI.COMM_WORLD, "scalar_wave.bp", [u_sol], engine="BP4") as vtx:
vtx.write(0.0)
