Hello. I am a new user of dolfinx and I am trying to convert the code in this link, from dolfin to dolfinx. The code in dolfinx is:
from mpi4py import MPI
from dolfinx import mesh, fem, plot
from dolfinx.fem import FunctionSpace, Function, Constant, Expression, dirichletbc, VectorFunctionSpace
import numpy as np
from petsc4py import PETSc
from petsc4py.PETSc import ScalarType
import ufl
from ufl import TrialFunction, TestFunction
import pyvista
import dolfinx
import gmsh
from dolfinx.io import gmshio
gmsh.initialize()
Sphere = gmsh.model.occ.addSphere(0, 0, 0, 1, 1)
gmsh.model.occ.synchronize()
Geometry_Dimension = 3
volumes = gmsh.model.getEntities(Geometry_Dimension)
Volume_Physical_Tag = 1
Volume_Physical = gmsh.model.addPhysicalGroup(volumes[0][0], [volumes[0][1]], Volume_Physical_Tag)
gmsh.model.mesh.field.add("Ball", 1)
gmsh.model.mesh.field.setNumber(1, "VIn", 0.1)
gmsh.model.mesh.field.setNumber(1, "XCenter", 0)
gmsh.model.mesh.field.setNumber(1, "YCenter", 0)
gmsh.model.mesh.field.setNumber(1, "ZCenter", 0)
gmsh.model.mesh.field.setNumber(1, "Radius", 1)
gmsh.model.mesh.field.add("Min", 2)
gmsh.model.mesh.field.setNumbers(2, "FieldsList", [1])
gmsh.model.mesh.field.setAsBackgroundMesh(2)
gmsh.model.mesh.generate(3)
mesh_comm = MPI.COMM_WORLD
model_rank = 0
domain, cell_tags, facet_tags = gmshio.model_to_mesh(gmsh.model, mesh_comm, model_rank, Geometry_Dimension)
# Define function spaces
V = FunctionSpace(domain, ("Nedelec 1st kind H(curl)", 1))
# Define test and trial functions
v = TestFunction(V)
u = TrialFunction(V)
# Define functions for boundary condiitons
dbdt = Constant(domain, (0.0, 0.0, 1.0))
zero = fem.Function(V)
zero_val = fem.Constant(domain, PETSc.ScalarType([0.0, 0.0, 0.0]))
zero_expr = fem.Expression(zero_val, V.element.interpolation_points())
zero.interpolate(zero_expr)
# Magnetic field (to be computed)
T = Function(V)
# Dirichlet boundary condition
tdim = domain.topology.dim
fdim = tdim - 1
domain.topology.create_connectivity(fdim, tdim)
boundary_facets = mesh.exterior_facet_indices(domain.topology)
boundary_dofs = fem.locate_dofs_topological(V, fdim, boundary_facets)
bc = dirichletbc(zero, boundary_dofs)
# Forms for the eddy-current equation
a = ufl.inner(ufl.curl(v), ufl.curl(u))*ufl.dx
L = -ufl.inner(v, dbdt)*ufl.dx
a_assembled = fem.petsc.assemble_matrix(fem.form(a), bcs = [bc])
a_assembled.assemble()
L_assembled = fem.petsc.assemble_vector(fem.form(L))
L_assembled.assemble()
# Create PETSc Krylov solver (from petsc4py)
ksp = PETSc.KSP()
ksp.create(PETSc.COMM_WORLD)
# Set the Krylov solver type and set tolerances
ksp.setType("cg")
ksp.setTolerances(rtol=1.0e-08, atol=1.0e-12, divtol=1.0e10, max_it=300)
# Get the preconditioner and set type (HYPRE AMS)
pc = ksp.getPC()
pc.setType("hypre")
pc.setHYPREType("ams")
# Build discrete gradient
P1 = FunctionSpace(domain, ("CG", 1))
G = dolfinx.cpp.fem.petsc.discrete_gradient(P1._cpp_object, V._cpp_object)
G.assemble()
# Attach discrete gradient to preconditioner
pc.setHYPREDiscreteGradient(G)
# Build constants basis for the Nedelec space
constants = [fem.Function(V) for i in range(3)]
for i, c in enumerate(constants):
direction = Constant(domain, [1.0 if i == j else 0.0 for j in range(3)])
c_expr = Expression(direction , V.element.interpolation_points())
c.interpolate(c_expr)
# Inform preconditioner of constants in the Nedelec space
cvecs = [constant.vector for constant in constants]
pc.setHYPRESetEdgeConstantVectors(cvecs[0], cvecs[1], cvecs[2])
# We are dealing with a zero conductivity problem (no mass term), so
# we need to tell the preconditioner
pc.setHYPRESetBetaPoissonMatrix(None)
# Set operator for the linear solver
ksp.setOperators(a_assembled)
# Set options prefix
ksp.setOptionsPrefix("eddy_")
# Turn on monitoring of residual
opts = PETSc.Options()
opts.setValue("-eddy_ksp_monitor_true_residual", None)
# Solve eddy currents equation (using potential T)
ksp.setFromOptions()
ksp.solve(L_assembled, T.vector)
# Test and trial functions for density equation
W = VectorFunctionSpace(domain, ("CG", 1))
v2 = TestFunction(W)
u2 = TrialFunction(W)
# Solve density equation (use conjugate gradient linear solver)
J = Function(W)
a2 = ufl.inner(v2, u2)*ufl.dx
L2 = ufl.inner(v2, ufl.curl(T))*ufl.dx
problem = fem.petsc.LinearProblem(a2, L2, petsc_options={"linear_solver": "cg"})
J = problem.solve()
However, by plotting J, the result on the boundary nodes are not match with the result from dolfin. The current density resulted from dolfin is shown in the image below.
The result from dolfinx is also shown in the image below.
Although the current density calculated for inner elements seems to be identical, but the result on the boundary is not correct in the code of dolfinx. I am wondering what the problem is and how I can fix my code. I would be so grateful if you could please help me in this regard.
Thank you in advance for your guidance.