Well, the question was about assembling before setting the nullspace, which is not required. However, to test the nullspace, the matrix needs to be assembled.
1 Like
Iām back at it @dokken. I am too ignorant to properly assert the test of the nullspace (I tried several variants, involving assemble_matrix() and a few others, I couldnāt get the proper syntax).
Nevertheless, here is my attempt without Newton solver, just like you asked for:
import numpy as np
from dolfinx import log, fem
from dolfinx.fem import (Constant, Function, functionspace, FunctionSpace, assemble_scalar,
dirichletbc, form, locate_dofs_geometrical, VectorFunctionSpace,
locate_dofs_topological, Expression, create_matrix, assemble_matrix)
from dolfinx.fem.petsc import NonlinearProblem, LinearProblem
from dolfinx.nls.petsc import NewtonSolver
from mpi4py import MPI
from petsc4py import PETSc
from petsc4py.PETSc import ScalarType
from ufl import (FiniteElement, Measure, MixedElement, SpatialCoordinate,
TestFunction, TrialFunction, as_tensor, dot, dx, grad, inner,
inv, split, derivative, FacetNormal, div)
import gmsh
from dolfinx.io import gmshio, VTXWriter
mesh, cell_markers, facet_markers = gmshio.read_from_msh("meshes/rectangles.msh", MPI.COMM_WORLD, gdim=2)
def resistance(meshfile):
the_current = 1.0
current_curves=(12,13)
voltages_curves=(12,13)
inj_current_curve, out_current_curve = current_curves
reading_voltage_curve_0, reading_voltage_curve_1 = voltages_curves
# Define FE function space
deg = 3
el = FiniteElement("CG", mesh.ufl_cell(), deg)
V = FunctionSpace(mesh, el)
v = TestFunction(V)
volt = TrialFunction(V)
rho = 1
sigma = 1.0 / rho
# Define the boundary conditions
left_facets = facet_markers.find(inj_current_curve)
right_facets = facet_markers.find(out_current_curve)
left_dofs = locate_dofs_topological(V, mesh.topology.dim-1, left_facets)
V_current_contact_out = 0
bcs = [dirichletbc(ScalarType(V_current_contact_out), left_dofs, V)]
x = SpatialCoordinate(mesh)
dx = Measure("dx", domain=mesh, subdomain_data=cell_markers)
ds = Measure("ds", domain=mesh, subdomain_data=facet_markers)
J = the_current / assemble_scalar(form(1 * ds(inj_current_curve, domain=mesh)))
# Weak form.
J_vector = -sigma * grad(volt)
a = dot(grad(v), J_vector)*dx
L = - v * J * ds(out_current_curve) - v * J * ds(inj_current_curve)
print(f''' ------- Pre-processing --------
Current density: {J}
Length of the side where current is injected: {assemble_scalar(form(1 * ds(inj_current_curve, domain=mesh)))}
Length of the side where current leaves the wire: {assemble_scalar(form(1 * ds(out_current_curve, domain=mesh)))}
This should correspond to the current injected: {assemble_scalar(form(J*ds(out_current_curve)))}
''')
A = fem.petsc.assemble_matrix(fem.form(a))
A.assemble()
b = fem.petsc.assemble_vector(fem.form(L))
b.ghostUpdate(addv=PETSc.InsertMode.ADD_VALUES, mode=PETSc.ScatterMode.REVERSE)
# Create vector that spans the null space
nullspace = PETSc.NullSpace().create(constant=True,comm=MPI.COMM_WORLD)
assert nullspace.test(A)
A.setNullSpace(nullspace)
nullspace.remove(b)
ksp = PETSc.KSP().create(mesh.comm)
ksp.setOperators(A)
opts = PETSc.Options()
opts["ksp_type"] = "preonly"
opts["pc_type"] = "lu"
opts["pc_factor_solver_type"] = "mumps"
opts["mat_mumps_icntl_24"] = 1 # Option to support solving a singular matrix
opts["mat_mumps_icntl_25"] = 0 # Option to support solving a singular matrix
ksp.setFromOptions()
volth = fem.Function(V)
ksp.solve(b, volth.vector)
volth.x.scatter_forward()
iterations = ksp.getIterationNumber()
# See: https://petsc.org/release/manualpages/KSP/KSPConvergedReason/
converged_reason = ksp.getConvergedReason()
if converged_reason > 0:
print(f"PETSc solver converged in {iterations=}")
else:
raise RuntimeError(f"PETSc solver did not converge with {converged_reason=}")
#problem = LinearProblem(a, L, petsc_options={"ksp_type": "preonly", "pc_type": "lu", "pc_factor_mat_solver_type": "mumps", "mat_mumps_icntl_24": 1, "mat_mumps_icntl_25" :0})
avg_voltage_curve_integral_0 = assemble_scalar(form(volt * ds(reading_voltage_curve_0, domain=mesh))) / assemble_scalar(form(1 * ds(reading_voltage_curve_0, domain=mesh)))
avg_voltage_curve_integral_1 = assemble_scalar(form(volt * ds(reading_voltage_curve_1, domain=mesh))) / assemble_scalar(form(1 * ds(reading_voltage_curve_1, domain=mesh)))
voltage_diff = avg_voltage_curve_integral_1 - avg_voltage_curve_integral_0
resistance = voltage_diff / assemble_scalar(form(J*ds(out_current_curve)))
n = FacetNormal(mesh)
down_current_flux = form(dot(J_vector, n)*ds(inj_current_curve))
I_out = assemble_scalar(down_current_flux)
up_current_flux = form(dot(J_vector, n)*ds(out_current_curve))
I_in = assemble_scalar(up_current_flux)
# Positive means outwards the curve/surface.
print(f'I_out: {I_out}')
print(f'I_in: {I_in}')
print(f'div(Je) should be 0: {assemble_scalar(form(div(J_vector) * dx(domain=mesh)))}')
return resistance
res = resistance("meshes/rectangles.msh")
print(f'resistance: {res}')
Hereās the output:
Info : Reading 'meshes/rectangles.msh'...
Info : 21 entities
Info : 1291 nodes
Info : 2353 elements
Info : Done reading 'meshes/rectangles.msh'
------- Pre-processing --------
Current density: 33.333333333333336
Length of the side where current is injected: 0.03
Length of the side where current leaves the wire: 0.09999999999999998
This should correspond to the current injected: 3.3333333333333335
PETSc solver converged in iterations=1
Loguru caught a signal: SIGSEGV
Stack trace:
24 0x5eb6ad3dcba5 _start + 37
23 0x7b05c8ff8e40 __libc_start_main + 128
22 0x7b05c8ff8d90 /lib/x86_64-linux-gnu/libc.so.6(+0x29d90) [0x7b05c8ff8d90]
21 0x5eb6ad3dccad Py_BytesMain + 45
20 0x5eb6ad40641e Py_RunMain + 702
19 0x5eb6ad413653 _PyRun_AnyFileObject + 67
18 0x5eb6ad413a08 _PyRun_SimpleFileObject + 424
17 0x5eb6ad414525 python3(+0x264525) [0x5eb6ad414525]
16 0x5eb6ad40e0bb python3(+0x25e0bb) [0x5eb6ad40e0bb]
15 0x5eb6ad4147d8 python3(+0x2647d8) [0x5eb6ad4147d8]
14 0x5eb6ad3e9cf6 PyEval_EvalCode + 134
13 0x5eb6ad3e9e56 python3(+0x239e56) [0x5eb6ad3e9e56]
12 0x5eb6ad2f8e0d _PyEval_EvalFrameDefault + 1725
11 0x5eb6ad31070c _PyFunction_Vectorcall + 124
10 0x5eb6ad2f8e0d _PyEval_EvalFrameDefault + 1725
9 0x5eb6ad31070c _PyFunction_Vectorcall + 124
8 0x5eb6ad2ff1f1 _PyEval_EvalFrameDefault + 27297
7 0x5eb6ad3065eb _PyObject_MakeTpCall + 603
6 0x5eb6ad30fe0e python3(+0x15fe0e) [0x5eb6ad30fe0e]
5 0x7b05bd377c5c /usr/local/dolfinx-real/lib/python3.10/dist-packages/dolfinx/cpp.cpython-310-x86_64-linux-gnu.so(+0xb1c5c) [0x7b05bd377c5c]
4 0x7b05bd3b2e0d /usr/local/dolfinx-real/lib/python3.10/dist-packages/dolfinx/cpp.cpython-310-x86_64-linux-gnu.so(+0xece0d) [0x7b05bd3b2e0d]
3 0x7b05bd62afbc /usr/local/dolfinx-real/lib/python3.10/dist-packages/dolfinx/cpp.cpython-310-x86_64-linux-gnu.so(+0x364fbc) [0x7b05bd62afbc]
2 0x7b05bd62c063 /usr/local/dolfinx-real/lib/python3.10/dist-packages/dolfinx/cpp.cpython-310-x86_64-linux-gnu.so(+0x366063) [0x7b05bd62c063]
1 0x7b05bd37e175 /usr/local/dolfinx-real/lib/python3.10/dist-packages/dolfinx/cpp.cpython-310-x86_64-linux-gnu.so(+0xb8175) [0x7b05bd37e175]
0 0x7b05c9011520 /lib/x86_64-linux-gnu/libc.so.6(+0x42520) [0x7b05c9011520]
2024-04-18 20:22:47.085 ( 0.841s) [main ] :0 FATL| Signal: SIGSEGV
Segmentation fault (core dumped)
A core dump, but at least the code converges within a single step, right? Unfortunately I donāt have the feedback of my print statements, to see if the values are OK.
By debugging by adding print statements, I figured out that this is the line that produces the core dump:
I_out = assemble_scalar(down_current_flux)
Edit: I could get those quantities:
avg_voltage_curve_integral_0: 0.08333333333333334
avg_voltage_curve_integral_1: 0.08333333333333336
voltage_diff: 1.3877787807814457e-17
resistance: 4.163336342344337e-18
It looks like it isnāt solved correctly, the correct resistance is near unity, so thereās a much bigger voltage difference than the one computed here.
Some more insights. With Paraview, I could examine the solution (after saving it with VTXwriter). It looks correct to me:
With bounds [-12.7152832665543, 15.403668747733219], this makes the resistance worth about 28 V /4 A, approximately 7 ohm, this is what I should get with FEniCSx.
Therefore, none of the 2 values in
avg_voltage_curve_integral_0: 0.08333333333333334
avg_voltage_curve_integral_1: 0.08333333333333336
are correct. And I canāt get a visual of the current density, as the code returns an error:
Traceback (most recent call last):
File "/root/shared/TE_simulation/electric_pb.py", line 195, in <module>
res = resistance("meshes/rectangles.msh")
File "/root/shared/TE_simulation/electric_pb.py", line 163, in resistance
Je.interpolate(Je_flux_calculator)
File "/usr/local/dolfinx-real/lib/python3.10/dist-packages/dolfinx/fem/function.py", line 397, in interpolate
_interpolate(u, cells)
File "/usr/lib/python3.10/functools.py", line 889, in wrapper
return dispatch(args[0].__class__)(*args, **kw)
File "/usr/local/dolfinx-real/lib/python3.10/dist-packages/dolfinx/fem/function.py", line 388, in _
self._cpp_object.interpolate(expr._cpp_object, cells)
RuntimeError: Cannot interpolate Expression with Argument
WARNING! There are options you set that were not used!
WARNING! could be spelling mistake, etc!
There are 3 unused database options. They are:
Option left: name:-mat_mumps_icntl_24 value: 1 source: code
Option left: name:-mat_mumps_icntl_25 value: 0 source: code
Option left: name:-pc_factor_solver_type value: mumps source: code
Soā¦ something is very strange. It looks like it correctly computed the solution (I suppose), but canāt evaluate anything on the boundaries?
Great @dokken, I could make some progress. It turns out that unlike when I use a Nonlinearsolver, here I must update J_vector, it isnāt updated automatically, hence the core dumps.
I fixed the code. It is now:
import numpy as np
from dolfinx import log, fem
from dolfinx.fem import (Constant, Function, functionspace, FunctionSpace, assemble_scalar,
dirichletbc, form, locate_dofs_geometrical, VectorFunctionSpace,
locate_dofs_topological, Expression, create_matrix, assemble_matrix)
from dolfinx.fem.petsc import NonlinearProblem, LinearProblem
from dolfinx.nls.petsc import NewtonSolver
from mpi4py import MPI
from petsc4py import PETSc
from petsc4py.PETSc import ScalarType
from ufl import (FiniteElement, Measure, MixedElement, SpatialCoordinate,
TestFunction, TrialFunction, as_tensor, dot, dx, grad, inner,
inv, split, derivative, FacetNormal, div)
import gmsh
from dolfinx.io import gmshio, VTXWriter
mesh, cell_markers, facet_markers = gmshio.read_from_msh("meshes/rectangles.msh", MPI.COMM_WORLD, gdim=2)
def resistance(meshfile):
the_current = 1.0
current_curves=(12,13)
voltages_curves=(12,13)
inj_current_curve, out_current_curve = current_curves
reading_voltage_curve_0, reading_voltage_curve_1 = voltages_curves
# Define FE function space
deg = 3
el = FiniteElement("CG", mesh.ufl_cell(), deg)
V = FunctionSpace(mesh, el)
v = TestFunction(V)
volt = TrialFunction(V)
rho = 1
sigma = 1.0 / rho
# Define the boundary conditions
left_facets = facet_markers.find(inj_current_curve)
right_facets = facet_markers.find(out_current_curve)
left_dofs = locate_dofs_topological(V, mesh.topology.dim-1, left_facets)
V_current_contact_out = 0
x = SpatialCoordinate(mesh)
dx = Measure("dx", domain=mesh, subdomain_data=cell_markers)
ds = Measure("ds", domain=mesh, subdomain_data=facet_markers)
J = the_current / assemble_scalar(form(1 * ds(inj_current_curve, domain=mesh)))
# Weak form.
J_vector = -sigma * grad(volt)
a = -dot(grad(v), J_vector)*dx
L = - v * J * ds(out_current_curve) + v * J * ds(inj_current_curve)
print(f''' ------- Pre-processing --------
Current density: {J}
Length of the side where current is injected: {assemble_scalar(form(1 * ds(inj_current_curve, domain=mesh)))}
Length of the side where current leaves the wire: {assemble_scalar(form(1 * ds(out_current_curve, domain=mesh)))}
This should correspond to the current injected: {assemble_scalar(form(J*ds(out_current_curve)))}
''')
A = fem.petsc.assemble_matrix(fem.form(a))
A.assemble()
b = fem.petsc.assemble_vector(fem.form(L))
b.ghostUpdate(addv=PETSc.InsertMode.ADD_VALUES, mode=PETSc.ScatterMode.REVERSE)
# Create vector that spans the null space
nullspace = PETSc.NullSpace().create(constant=True,comm=MPI.COMM_WORLD)
assert nullspace.test(A)
A.setNullSpace(nullspace)
nullspace.remove(b)
ksp = PETSc.KSP().create(mesh.comm)
ksp.setOperators(A)
opts = PETSc.Options()
opts["ksp_type"] = "preonly"
opts["pc_type"] = "lu"
opts["pc_factor_solver_type"] = "mumps"
opts["mat_mumps_icntl_24"] = 1 # Option to support solving a singular matrix
opts["mat_mumps_icntl_25"] = 0 # Option to support solving a singular matrix
ksp.setFromOptions()
volth = fem.Function(V)
ksp.solve(b, volth.vector)
volth.x.scatter_forward()
iterations = ksp.getIterationNumber()
# See: https://petsc.org/release/manualpages/KSP/KSPConvergedReason/
converged_reason = ksp.getConvergedReason()
if converged_reason > 0:
print(f"PETSc solver converged in {iterations=}")
else:
raise RuntimeError(f"PETSc solver did not converge with {converged_reason=}")
with VTXWriter(MPI.COMM_WORLD, "results/voltage_elec.bp", [volth], engine="BP4") as vtx:
vtx.write(0.0)
# Update J_vector
J_vector = -sigma * grad(volth)
avg_voltage_curve_integral_0 = assemble_scalar(form(volth * ds(reading_voltage_curve_0, domain=mesh))) / assemble_scalar(form(1 * ds(reading_voltage_curve_0, domain=mesh)))
avg_voltage_curve_integral_1 = assemble_scalar(form(volth * ds(reading_voltage_curve_1, domain=mesh))) / assemble_scalar(form(1 * ds(reading_voltage_curve_1, domain=mesh)))
voltage_diff = avg_voltage_curve_integral_1 - avg_voltage_curve_integral_0
resistance = voltage_diff / assemble_scalar(form(J*ds(out_current_curve)))
print(f'avg_voltage_curve_integral_0: {avg_voltage_curve_integral_0}')
print(f'avg_voltage_curve_integral_1: {avg_voltage_curve_integral_1}')
print(f'voltage_diff: {voltage_diff}')
print(f'resistance: {resistance}')
print(f'div(Je) should be 0: {assemble_scalar(form(div(J_vector) * dx(domain=mesh)))}')
n = FacetNormal(mesh)
down_current_flux = form(dot(J_vector, n)*ds(inj_current_curve))
I_out = assemble_scalar(down_current_flux)
up_current_flux = form(dot(J_vector, n)*ds(out_current_curve))
I_in = assemble_scalar(up_current_flux)
# Positive means outwards the curve/surface.
print(f'I_out: {I_out}')
print(f'I_in: {I_in}')
return resistance
res = resistance("meshes/rectangles.msh")
print(f'resistance: {res}')
The output is:
Info : Reading 'meshes/rectangles.msh'...
Info : 21 entities
Info : 1291 nodes
Info : 2353 elements
Info : Done reading 'meshes/rectangles.msh'
------- Pre-processing --------
Current density: 33.333333333333336
Length of the side where current is injected: 0.03
Length of the side where current leaves the wire: 0.09999999999999998
This should correspond to the current injected: 3.3333333333333335
PETSc solver converged in iterations=1
avg_voltage_curve_integral_0: 15.403275750479303
avg_voltage_curve_integral_1: -12.714402930356114
voltage_diff: -28.117678680835418
resistance: -8.435303604250626
div(Je) should be 0: 1.1140697211743482
I_out: -1.0014547932721676
I_in: 3.3286405396011376
resistance: -8.435303604250626
WARNING! There are options you set that were not used!
WARNING! could be spelling mistake, etc!
There are 3 unused database options. They are:
Option left: name:-mat_mumps_icntl_24 value: 1 source: code
Option left: name:-mat_mumps_icntl_25 value: 0 source: code
Option left: name:-pc_factor_solver_type value: mumps source: code
Which points at 2 problems (maybe related):
- For some reason, the solver options to take into account the nullspace are ignored.
- The solution computed is wrong. I am solving div(J)=0 and I get a value of 1.11. The currents in and out are also wrong, the in current is slightly off, the out current is quite off the mark.
Thatās the same kind of problems I reported in post #9.
Any idea what could be wrong? I suspect the nullspace isnāt taken well into account. I donāt see anything else at the moment.
Edit: I am trying ārandomā stuff. I could get to use mumps (I believe, though it doesnāt change any result so far) by adding the line ksp.getPC().setFactorSolverType("MUMPS"), notice the capital letters (Iād get the warning that 3 options were not used, as above, otherwise). With these capital letters, I only have the remaining 2 warnings to get rid ofā¦ if anyone reads me, feel free to help. Iām desperately trying to get this to work properly.
Edit2: Trying the syntax:
ksp.getPC().getFactorMatrix().setMumpsIcntl(icntl=24, ival=1)
ksp.getPC().getFactorMatrix().setMumpsIcntl(icntl=25, ival=0)
yields:
File "/root/shared/TE_simulation/electric_pb.py", line 187, in <module>
res = resistance("meshes/rectangles.msh")
File "/root/shared/TE_simulation/electric_pb.py", line 106, in resistance
ksp.getPC().getFactorMatrix().setMumpsIcntl(icntl=24, ival=1)
File "petsc4py/PETSc/PC.pyx", line 1433, in petsc4py.PETSc.PC.getFactorMatrix
File "petsc4py/PETSc/PETSc.pyx", line 79, in petsc4py.PETSc.CHKERR
petsc4py.PETSc.Error: error code 56
[0] PCFactorGetMatrix() at /usr/local/petsc/src/ksp/pc/interface/precon.c:1456
[0] No support for this operation for this object type
[0] No method getfactoredmatrix for PC of type (null)
At least the statements are recognizedā¦ though they lead to an errorā¦
I am still totally stuck on this. If someone knows how I can use mumps with the null space options as specified, let me know. I tried at least 3 different syntaxes, and so far it seems only 1 works for mumps but I havenāt been able to set the null space options.
Iām not sure people are still reading me.
I just wanted to mention that chatGPT helped me an iota, and pointed out that my opts["pc_factor_solver_type"] = "mumps" should in fact have been opts["pc_factor_mat_solver_type"] = "mumps"
So now I know mumps is used, with the options I have set. However this doesnāt change the wrong result, I get exactly the same wrong result.
Hereās some ādebugā info from mumps:
Entering DMUMPS 5.6.1 from C interface with JOB, N, NNZ = 1 10891 180811
executing #MPI = 1, without OMP
=================================================
MUMPS compiled with option -Dmetis
MUMPS compiled with option -Dpord
MUMPS compiled with option -Dptscotch
MUMPS compiled with option -Dscotch
=================================================
L U Solver for unsymmetric matrices
Type of parallelism: Working host
****** ANALYSIS STEP ********
Processing a graph of size: 10891
Entering ordering phase with ...
N NNZ LIW INFO(1)
10891 180811 372514 0
Matrix entries: IRN() ICN()
1 1 1 2 1 3
1 4 1 5 1 6
1 7 1 8 1 9
1 10
... Structural symmetry (in percent)= 100
Average density of rows/columns = 16
... No column permutation
Ordering based on METIS
ELAPSED TIME SPENT IN METIS reordering = 0.0333
SYMBOLIC based on column counts
ELAPSED TIME IN symbolic factorization = 321.2459
NFSIZ(.) = 0 0 0 0 31 0 0 0 0 0
FILS (.) = 12 24 2 3 27 0 6 1 8 7
FRERE(.) = 10892 10892 10892 10892 -54 10892 10892 10892 10892 10892
Leaving analysis phase with ...
INFOG(1) = 0
INFOG(2) = 0
-- (20) Number of entries in factors (estim.) = 845915
-- (3) Real space for factors (estimated) = 845915
-- (4) Integer space for factors (estimated) = 73224
-- (5) Maximum frontal size (estimated) = 118
-- (6) Number of nodes in the tree = 688
-- (32) Type of analysis effectively used = 1
-- (7) Ordering option effectively used = 5
ICNTL (6) Maximum transversal option = 0
ICNTL (7) Pivot order option = 7
ICNTL(13) Parallelism/splitting of root node = 1
ICNTL(14) Percentage of memory relaxation = 20
ICNTL(15) Analysis by block effectively used = 0
ICNTL(18) Distributed input matrix (on if >0) = 0
ICNTL(58) Symbolic factorization option = 2
Number of level 2 nodes = 0
Number of split nodes = 0
RINFOG(1) Operations during elimination (estim)= 3.975D+07
MEMORY ESTIMATIONS ...
Estimations with standard Full-Rank (FR) factorization:
Total space in MBytes, IC factorization (INFOG(17)): 11
Total space in MBytes, OOC factorization (INFOG(27)): 4
Elapsed time in analysis driver= 0.0397
Entering DMUMPS 5.6.1 from C interface with JOB, N, NNZ = 2 10891 180811
executing #MPI = 1, without OMP
****** FACTORIZATION STEP ********
GLOBAL STATISTICS PRIOR NUMERICAL FACTORIZATION ...
Number of working processes = 1
ICNTL(22) Out-of-core option = 0
ICNTL(35) BLR activation (eff. choice) = 0
ICNTL(37) BLR CB compression (eff. choice) = 0
ICNTL(49) Compact workarray S (end facto.) = 0
ICNTL(14) Memory relaxation = 20
INFOG(3) Real space for factors (estimated)= 845915
INFOG(4) Integer space for factors (estim.)= 73224
Maximum frontal size (estimated) = 118
Number of nodes in the tree = 688
ICNTL(23) Memory allowed (value on host) = 0
Sum over all procs = 0
Memory provided by user, sum of LWK_USER = 0
Effective threshold for pivoting, CNTL(1) = 0.1000D-01
Max difference from 1 after scaling the entries for ONE-NORM (option 7/8) = 0.40D-01
CNTL(3) for null pivot row detection = 0.0000D+00
ICNTL(24) null pivot rows detection = 1
...Zero pivot detection threshold = 4.5257D-15
...Infinite fixation
Effective size of S (based on INFO(39))= 1024969
Elapsed time to reformat/distribute matrix = 0.0029
** Memory allocated, total in Mbytes (INFOG(19)): 11
** Memory effectively used, total in Mbytes (INFOG(22)): 10
Elapsed time for factorization = 0.0110
Leaving factorization with ...
RINFOG (2) Operations in node assembly = 6.493D+05
------ (3) Operations in node elimination = 3.975D+07
ICNTL (8) Scaling effectively used = 7
INFOG (9) Real space for factors = 845915
INFOG (10) Integer space for factors = 73224
INFOG (11) Maximum front size = 118
INFOG (29) Number of entries in factors = 845915
INFOG (12) Number of off diagonal pivots = 0
INFOG (13) Number of delayed pivots = 0
Nb of null pivots detected by ICNTL(24) = 0
INFOG (28) Estimated deficiency = 0
INFOG (14) Number of memory compress = 0
Elapsed time in factorization driver = 0.0172
Entering DMUMPS 5.6.1 from C interface with JOB, N, NNZ = 3 10891 180811
executing #MPI = 1, without OMP
****** SOLVE & CHECK STEP ********
GLOBAL STATISTICS PRIOR SOLVE PHASE ...........
Number of right-hand-sides = 1
Blocking factor for multiple rhs = 1
ICNTL (9) = 1
--- (10) = 0
--- (11) = 0
--- (20) = 0
--- (21) = 0
--- (30) = 0
--- (35) = 0
Vector solution for column 1
RHS
1.432576D+01 1.408039D+01 1.399049D+01 1.405549D+01 1.401528D+01
1.423322D+01 1.408331D+01 1.425793D+01 1.414823D+01 1.413229D+01
** Space in MBYTES used for solve : 9
Leaving solve with ...
Time to build/scatter RHS = 0.000108
Time in solution step (fwd/bwd) = 0.001701
.. Time in forward (fwd) step = 0.000907
.. Time in backward (bwd) step = 0.000793
Time to gather solution(cent.sol)= 0.000015
Time to copy/scale dist. solution= 0.000000
Elapsed time in solve driver= 0.0022
Edit: You can consider this thread as finished. I will post my question in a new thread, since the original question has been answered, I can get a convergence in 1 step using either 2 Dirichlet b.c. or a combination of 1 Dirichlet and 1 Neumann b.c.