Magnetostatics 3D - wrong 3D magnetic vector potential

General
5

Thanks for the quick answer and the link. I made the change.

I could not get ADIOS2 to compile correctly on my ubuntu 22.04LTS machine to get VTXWriter working. It seems to be an issue of debian systems.

from mpi4py import MPI
from dolfinx.io import gmshio 
from dolfinx.io import XDMFFile
from dolfinx.fem import (Expression, Function, FunctionSpace, VectorFunctionSpace)
from dolfinx.fem import (dirichletbc, locate_dofs_topological)
from dolfinx import *
from dolfinx import io

import ufl
from ufl import as_vector, dot, dx, grad, curl, inner

from petsc4py.PETSc import ScalarType
import numpy as np

rank = MPI.COMM_WORLD.rank

import gmsh
import sys
import math

gmsh.initialize()

gmsh.model.add("magnetostatics3D")

gmsh.option.setColor("Mesh.Color.Lines", 0, 0, 0)

model_rank = 0
mesh_comm = MPI.COMM_WORLD
gdim = 3

Ri = 0.1
th = 0.05
Ra = 0.3
la = 0.5

lc = th / 10
nl = 60

if mesh_comm.rank == model_rank:
	cascade = gmsh.model.occ
	meshing = gmsh.model.mesh

	cascade.addPoint(Ri, 0, -th/2, lc, 1)
	cascade.addPoint(Ri+th, 0, -th/2, lc, 2)
	cascade.addPoint(Ri+th, 0, th/2, lc, 3)
	cascade.addPoint(Ri, 0, th/2, lc, 4)
	cascade.addLine(1, 2, 1)
	cascade.addLine(2, 3, 2)
	cascade.addLine(3, 4, 3)
	cascade.addLine(4, 1, 4)
	cascade.addCurveLoop([1, 2, 3, 4], 1)
	cascade.addPlaneSurface([1], 1)
	
	cascade.addPoint(0., 0., -la/2, lc, 5)
	cascade.addPoint(0., 0., la/2, lc, 6)
	cascade.addLine(5, 6, 5)
	
	cascade.synchronize()

	ls_1 = cascade.revolve([(2, 1)], 0, 0, 0, 0, 0, 1, -math.pi, [nl], recombine = False)
	ls_2 = cascade.revolve([(2, 1)], 0, 0, 0, 0, 0, 1, math.pi, [nl], recombine = False)
		
	cascade.fragment(gmsh.model.occ.getEntities(3), [])
	
	cascade.addSphere(0, 0, 0, Ra, 3)
	# Cut out the coil from the sphere keeping a unique interface for both
	cascade.fragment([(3, 3)], [(3, 1) ,(3, 2), (1, 5)])
	
	cascade.removeAllDuplicates
	meshing.removeDuplicateNodes
	meshing.removeDuplicateElements
	
	cascade.synchronize()
	
	up, down = gmsh.model.getAdjacencies(gdim, 3)

	gmsh.model.addPhysicalGroup(gdim, [1, 2], 1, name="Coil")
	gmsh.model.setColor(ls_1, 255, 0, 0)
	gmsh.model.setColor(ls_2, 255, 0, 0)
	gmsh.model.addPhysicalGroup(gdim, [3], 2, name="Air")
	gmsh.model.setColor([(gdim, 3)], 0, 0, 255)
	gmsh.model.addPhysicalGroup(gdim-1, [down[0]], 3, name="Boundary")
	gmsh.model.setColor([(gdim-1, down[0])], 0, 0, 255)
	gmsh.model.addPhysicalGroup(1, [5], 4, name="Axis")

	meshing.generate(gdim)

	gmsh.write("magnetostatics3D.msh")

# Read in by dolfinx
mesh, cell_tags, facet_tags = gmshio.model_to_mesh(gmsh.model, mesh_comm, model_rank, gdim=gdim)
gmsh.finalize()

# Output DOLFINx meshes to file
with XDMFFile(MPI.COMM_WORLD, "meshCheck.xdmf", "w") as xdmf:
    xdmf.write_mesh(mesh)
    xdmf.write_meshtags(facet_tags)
    xdmf.write_meshtags(cell_tags)

# Function space: edge elements
V0 = FunctionSpace(mesh, ("N1curl", 1))

# Create facet to cell connectivity required to determine boundary facets
dofs_facets = locate_dofs_topological(V0, mesh.topology.dim-1, facet_tags.find(3))

# Working boundary condition
u0 = Function(V0)
u0.x.set(0.)
bc = dirichletbc(u0, dofs_facets)

# Trial and test functions
u = ufl.TrialFunction(V0)
v = ufl.TestFunction(V0)

##Define Current Density and Mu0
V1 = VectorFunctionSpace(mesh, ("CG", 1))
J = Function(V1, dtype=ScalarType, name='Current_density')
cells1 = cell_tags.find(1)

J0 = 1e8
def Je(x, J0):
    return np.vstack([-J0*np.sin(np.arctan2(x[1],x[0]))*np.ones(x.shape[1]), J0*np.cos(np.arctan2(x[1],x[0]))*np.ones(x.shape[1]), np.zeros(x.shape[1])])

J.interpolate(lambda x: Je(x, J0), cells1)
J.x.scatter_forward()

mu = 4*np.pi*1e-7

a = (1 / mu) * inner(curl(u), curl(v)) * dx
l = inner(J, v) * dx

A = Function(V0)
problem = fem.petsc.LinearProblem(a, l, u = A, bcs = [bc], petsc_options = {"ksp_type": "gmres", "ksp_rtol":1e-8, "ksp_atol":1e-12, "ksp_max_it": 1000, "pc_type": "none"})
problem.solve()

V2 = VectorFunctionSpace(mesh, ("DG", 1))
B = Function(V2)
B_expr = Expression(ufl.curl(A), V2.element.interpolation_points())
B.interpolate(B_expr)

with io.XDMFFile(mesh.comm, "J.xdmf", "w") as xdmf:
    xdmf.write_mesh(mesh)
    xdmf.write_meshtags(cell_tags) 
    xdmf.write_function(J)

with io.XDMFFile(mesh.comm, "A.xdmf", "w") as xdmf:
    xdmf.write_mesh(mesh)
    xdmf.write_meshtags(cell_tags) 
    xdmf.write_function(A)
    
with io.XDMFFile(mesh.comm, "B.xdmf", "w") as xdmf:
    xdmf.write_mesh(mesh)
    xdmf.write_meshtags(cell_tags) 
    xdmf.write_function(B)

#with dolfinx.io.VTXWriter(mesh.comm, "B.bp", [B]) as vtx:
#    vtx.write(0.0)

Best,

Frederic