Hello FEniCSxers,
Tl;dr: jump at the last code and try to make it converge to correct values.
I am tackling a very simple PDE: \nabla \cdot \vec J=0 where \vec J=-\sigma \nabla V (Ohm’s law), which means that the temporal change in electric charge in any region vanishes, i.e. no creation/destruction of electric charge and any electric charge that enters any region must leave it, thereby creating an electric current.
It is a linear problem in 2D. The mesh can be generated with gmsh with this file.geo:
SetFactory("OpenCASCADE");
Mesh.CharacteristicLengthMax = 0.010;
Mesh.CharacteristicLengthMin = 0.001;
Rectangle(1) = {0, 0, 0, 0.03, 0.3, 0};
Rectangle(2) = {0, 0.2, 0, 0.9, 0.1, 0};
BooleanFragments{ Surface{2}; Surface{1}; Delete; }{ }
BooleanFragments{ Surface{1}; Surface{3}; Delete; }{ }
Physical Surface("material", 11) = {2, 1, 3};
Physical Curve("bottom_left", 12) = {10};
Physical Curve("top_right", 13) = {6};
using the GUI for example, and clicking on Save mesh.
I get a convergence in 1 step, as it should, when I use at least one Dirichlet boundary conditions. The full code is then:
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(inj_current_curve) - v * J * ds(out_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)))}
''')
# Solve the PDE.
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)
default_problem = LinearProblem(a, L, bcs=bcs, petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
volth = default_problem.solve()
volth.x.scatter_forward()
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 = abs(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_in = assemble_scalar(down_current_flux)
up_current_flux = form(dot(J_vector, n)*ds(out_current_curve))
I_out = assemble_scalar(up_current_flux)
# Positive means outwards the curve/surface.
print(f'I_in: {I_in}')
print(f'I_out: {I_out}')
# Current density.
J = functionspace(mesh, ("DG", deg-1, (mesh.geometry.dim,)))
Je = Function(J)
Je_flux_calculator = Expression(J_vector, J.element.interpolation_points())
Je.interpolate(Je_flux_calculator)
with VTXWriter(MPI.COMM_WORLD, "results/current_density_electric.bp", [Je], engine="BP4") as vtx:
vtx.write(0.0)
return resistance
res = resistance("meshes/rectangles.msh")
print(f'resistance: {res}')
The code returns:
Info : Reading 'meshes/rectangles.msh'...
Info : 21 entities
Info : 12931 nodes
Info : 25102 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.333333333333333
avg_voltage_curve_integral_0: 0.0
avg_voltage_curve_integral_1: -53.89244623146631
voltage_diff: -53.89244623146631
resistance: 16.167733869439893
div(Je) should be 0: 7.39213122723939e-09
I_in: -3.3333333331551662
I_out: 3.33333333334904
resistance: 16.167733869439893
which seem to be correct values. It is extremely important that the input electric current is equal to the output one, and that the divergence of the current density vanishes, after all, this is the equation I am solving, so if this condition is violated, it would mean that the equation has not been solved correctly. Note that I also compute the current density and save it for Paraview:
You can probably barely see it, but current goes in from the bottom left (red part, higher magnitude), and escapes the material on the right side.
All seems to be fine to me, so far.
But it turns out that my real problem is more complicated than this one, and I cannot apply a Dirichlet b.c. on a whole surface (at best, I could apply it at a node, if I had to). It turns out that Neumann b.c. are better suited for my problem, because they would correspond to what really happens to a sample.
So I should be able to solve this same problem using 2 Neumann b.c. (pure Neumann b.c.s). The problem is that there’s no unique solution in that case, voltage being defined up to an arbitrary constant, I must either specify the null space to the solver (that’s the way I wish to do things), or specify a voltage value at a node (I didn’t try yet, but I’d rather not do this, I prefer the solver do its thing).
Here’s the equivalent problem and code using pure Neumann b.c.s:
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(inj_current_curve) - v * J * ds(out_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)))}
''')
# Solve the PDE.
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)
# For direct solvers.
A.setNullSpace(nullspace)
# For iterative solvers.
#A.setNearNullSpace(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_mat_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
opts["mat_mumps_icntl_4"] = 3 # level of printing, from 0 to 4.
opts = PETSc.Options()
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 = abs(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_in = assemble_scalar(down_current_flux)
up_current_flux = form(dot(J_vector, n)*ds(out_current_curve))
I_out = assemble_scalar(up_current_flux)
# Positive means outwards the curve/surface.
print(f'I_in: {I_in}')
print(f'I_out: {I_out}')
# Current density.
J = functionspace(mesh, ("DG", deg-1, (mesh.geometry.dim,)))
Je = Function(J)
Je_flux_calculator = Expression(J_vector, J.element.interpolation_points())
Je.interpolate(Je_flux_calculator)
with VTXWriter(MPI.COMM_WORLD, "results/current_density_electric.bp", [Je], engine="BP4") as vtx:
vtx.write(0.0)
return resistance
res = resistance("meshes/rectangles.msh")
print(f'resistance: {res}')
The output is:
Info : Reading 'meshes/rectangles.msh'...
Info : 21 entities
Info : 12931 nodes
Info : 25102 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.333333333333333
Entering DMUMPS 5.6.1 from C interface with JOB, N, NNZ = 1 113965 1922953
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: 113965
Entering ordering phase with ...
N NNZ LIW INFO(1)
113965 1922953 3959872 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.5951
SYMBOLIC based on column counts
ELAPSED TIME IN symbolic factorization = 3195.5498
NFSIZ(.) = 0 0 0 0 0 0 0 0 0 0
FILS (.) = 25 41 2 3 4 10 0 1 8 7
FRERE(.) = 113966 113966 113966 113966 113966 113966 113966 113966 113966 113966
Leaving analysis phase with ...
INFOG(1) = 0
INFOG(2) = 0
-- (20) Number of entries in factors (estim.) = 11733651
-- (3) Real space for factors (estimated) = 11733651
-- (4) Integer space for factors (estimated) = 855768
-- (5) Maximum frontal size (estimated) = 363
-- (6) Number of nodes in the tree = 7148
-- (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)= 1.062D+09
MEMORY ESTIMATIONS ...
Estimations with standard Full-Rank (FR) factorization:
Total space in MBytes, IC factorization (INFOG(17)): 146
Total space in MBytes, OOC factorization (INFOG(27)): 38
Elapsed time in analysis driver= 0.6515
Entering DMUMPS 5.6.1 from C interface with JOB, N, NNZ = 2 113965 1922953
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)= 11733651
INFOG(4) Integer space for factors (estim.)= 855768
Maximum frontal size (estimated) = 363
Number of nodes in the tree = 7148
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.42D-01
CNTL(3) for null pivot row detection = 0.0000D+00
ICNTL(24) null pivot rows detection = 1
...Zero pivot detection threshold = 8.4847D-15
...Infinite fixation
Effective size of S (based on INFO(39))= 14186737
Elapsed time to reformat/distribute matrix = 0.0288
** Memory allocated, total in Mbytes (INFOG(19)): 146
** Memory effectively used, total in Mbytes (INFOG(22)): 127
Elapsed time for factorization = 0.1384
Leaving factorization with ...
RINFOG (2) Operations in node assembly = 1.240D+07
------ (3) Operations in node elimination = 1.062D+09
ICNTL (8) Scaling effectively used = 7
INFOG (9) Real space for factors = 11733651
INFOG (10) Integer space for factors = 855768
INFOG (11) Maximum front size = 363
INFOG (29) Number of entries in factors = 11733651
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.1974
Entering DMUMPS 5.6.1 from C interface with JOB, N, NNZ = 3 113965 1922953
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.644367D+01 -1.637027D+01 -1.644373D+01 -1.639057D+01 -1.642342D+01
-1.644370D+01 -1.644372D+01 -1.642339D+01 -1.639056D+01 -1.641923D+01
** Space in MBYTES used for solve : 126
Leaving solve with ...
Time to build/scatter RHS = 0.000630
Time in solution step (fwd/bwd) = 0.020997
.. Time in forward (fwd) step = 0.009877
.. Time in backward (bwd) step = 0.011117
Time to gather solution(cent.sol)= 0.000238
Time to copy/scale dist. solution= 0.000000
Elapsed time in solve driver= 0.0228
PETSc solver converged in iterations=1
avg_voltage_curve_integral_0: 15.390689284011646
avg_voltage_curve_integral_1: -12.557777313921909
voltage_diff: -27.948466597933553
resistance: 8.384539979380067
div(Je) should be 0: 1.1400901191042545
I_in: -1.00044703634129
I_out: 3.3318240875566643
Entering DMUMPS 5.6.1 from C interface with JOB = -2
executing #MPI = 1, without OMP
resistance: 8.384539979380067
Not only is \nabla \cdot \vec J \neq 0 (the PDE is not satisfied), but the current in is wrong (curiously the current out is correct).
The current density looks like:
which doesn’t make sense, as it is increasing along the length of the material instead of remaining roughly constant.
I haven’t been able to solve this simple problem with 2 Neumann b.c. so far. Any idea is welcome.