Hello everyone,
I tried to solve the problem of electric field distribution in a two-dimensional axisymmetric coordinate system (r-z), but I found that the output value of electric field strength is always “inf”.
Here is an MWE. The problem is a calculation domain with r=2mm and z=5mm, where the potential at the position of z=5mm is 10V and the potential at the position of z=0mm is 0.
import numpy
from dolfinx import mesh, fem
from dolfinx.fem.petsc import LinearProblem
import ufl
from ufl import (dot, dx, grad)
from petsc4py import PETSc
from mpi4py import MPI
domain = mesh.create_rectangle(MPI.COMM_WORLD, ([0,0],[2e-3,5e-3]), [20,50], mesh.CellType.triangle)
# Defining FunctionSpace
## Voltage
V = fem.FunctionSpace(domain, ("Lagrange", 1))
## Electric intensity
element_E = ufl.VectorElement("Lagrange", domain.ufl_cell(), 1)
E = fem.FunctionSpace(domain, element_E)
# Create boundary condition
## anode_marker
def anode_marker(x):
return numpy.isclose(x[1], 5e-3)
## cathode_marker
def cathode_marker(x):
return numpy.isclose(x[1], 0)
anode_boundary_entity = mesh.locate_entities_boundary(domain, 1, anode_marker)
cathode_boundary_entity = mesh.locate_entities_boundary(domain, 1, cathode_marker)
anode_boundary_dof = fem.locate_dofs_topological(V, 1, anode_boundary_entity)
bc1_voltage = fem.dirichletbc(PETSc.ScalarType(10), anode_boundary_dof, V)
cathode_boundary_dof = fem.locate_dofs_topological(V, 1, cathode_boundary_entity)
bc2_voltage = fem.dirichletbc(PETSc.ScalarType(0), cathode_boundary_dof, V)
# Define trailfunction and testfunction
## Voltage
u_voltage = ufl.TrialFunction(V)
v_voltage = ufl.TestFunction(V)
## electric intensity
u_E = ufl.TrialFunction(E)
v_E = ufl.TestFunction(E)
# Compute voltage
x = ufl.SpatialCoordinate(domain)
a = dot(grad(u_voltage), grad(v_voltage)) * x[0] * dx
f = fem.Constant(domain, PETSc.ScalarType(0))
L = f * v_voltage * x[0] * dx
problem = LinearProblem(a, L, bcs=[bc1_voltage, bc2_voltage], petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
u_voltage = problem.solve()
# Compute electric intensity
a = dot(u_E, v_E) * x[0] * dx
L = -dot(grad(u_voltage), v_E) * x[0] * dx
problem = LinearProblem(a, L, petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
u_E = problem.solve()
print(u_voltage.x.array)
print(u_E.x.array)
The result of running the code is as follows：
[ 0. 0. 0.2 ... 9.8 10. 10. ]
[inf inf inf ... inf inf inf]
It is obvious that the electric field strength generated by a uniform electric field with a potential difference of 10V and a distance of 5mm should be 2000V/m.
I have made some modifications that seem to avoid outputting “inf”.
If the number of grids is reduced, for example:
domain = mesh.create_rectangle(MPI.COMM_WORLD, ([0,0],[2e-3,5e-3]), [2,5], mesh.CellType.triangle)
The result of running the code is as follows：
[ 0. 0. 2. 2. 0. 4. 2. 4. 6. 4. 6. 8. 6. 8. 10. 8. 10. 10.]
[ 2.09704918e-13 -2.00000000e+03 -3.21492470e-13 -2.00000000e+03
4.56580570e-13 -2.00000000e+03 1.38844405e-13 -2.00000000e+03
-4.75498750e-13 -2.00000000e+03 3.45393962e-13 -2.00000000e+03
-7.44081027e-13 -2.00000000e+03 -2.87688412e-13 -2.00000000e+03
-4.77103187e-13 -2.00000000e+03 4.57738155e-13 -2.00000000e+03
-1.01184923e-12 -2.00000000e+03 -1.90857790e-12 -2.00000000e+03
-9.67609523e-13 -2.00000000e+03 -1.53526676e-12 -2.00000000e+03
-6.75744090e-13 -2.00000000e+03 -2.50548673e-12 -2.00000000e+03
1.44144609e-13 -2.00000000e+03 1.10859876e-12 -2.00000000e+03]
If I replace the coordinate system with a two-dimensional Cartesian coordinate system (x-y) — actually the only difference in the code is that x[0] is missing from the variational equation.
The result of running the code is as follows：
[ 0. 0. 0.2 ... 9.8 10. 10. ]
[-2.04011867e-13 -2.00000000e+03 8.99103355e-14 ... -2.00000000e+03
-1.17962349e-11 -2.00000000e+03]
These two appear to be correct，but what’s going on with the e-13 magnitude?
Therefore, there must be a problem in this, can you give me some instructions?