Hi all,
I am working on the curl curl problem, and trying to calculate the magnetic field around a torus filled with current in 3D. I am using the 0.6 version of dolfinx and used theses lines ot create my conda env
conda create -n vfenicsx -c conda-forge fenics-dolfinx
conda install -c conda-forge gmsh python-gmsh
I had some specific problem with N1curl even if in most of the case the magnetic field looked ok, so i try with N2curl. And then instead of having 1/r² field like it was supposed to be and what i had with N1curl i have something that decreases way faster.
my equation was ((1 / mu) * inner(rot(u), rot(v)) + 1e-6 * inner(u,v))* dx = inner(J, v) * dx
but i had to add to the left hand side the gauge “1e-6 * inner(u,v)*dx” for N2curl, it doesn’t change anything in N1curl though. Maybe it has something to do with it?
Here is a MWE (i did my best to make it short, sry, the code takes 10 second for N2curl and 10 to 20 times more for N1curl)
If anyone understand why these two look so different and how to correct the wrong one i"d appreciate your help. Thanks for reading me
from dolfinx.io import XDMFFile, gmshio
from mpi4py import MPI
import gmsh
import numpy as np
from mpi4py import MPI
import dolfinx
from dolfinx import io
from dolfinx import fem, mesh, plot
from dolfinx.fem import (Expression, VectorFunctionSpace, dirichletbc, Function, FunctionSpace, locate_dofs_geometrical,locate_dofs_topological)
from ufl import inner
from dolfinx.fem.petsc import LinearProblem
from ufl import TestFunction, TrialFunction, dx, inner, rot
### variable ###
refinement=5 # refinement of the meshing
space="N1curl" # in what space do we solve the equation
order=1 # what order of polynome
order_interpolation=1 # What order of interpolation to export the solution (in DG)
### Mesh Creation ###
mesh_name=f"{format(refinement, '.0e')}"
mesh_path=f"Mesh/{mesh_name}.xdmf"
gmsh.initialize()
gmsh.option.setNumber("General.Terminal", 0)
model = gmsh.model()
sphere_dim_tags = model.occ.addSphere(0, 0, 0, 100)
torus_dim_tag = model.occ.addTorus(0, 0, 0, 12.5, 1.5)
model_dim_tags = model.occ.fragment([(3, sphere_dim_tags)], [(3, torus_dim_tag)])
model.occ.synchronize()
# Add physical tag 1 for exterior surfaces (for bc)
boundary = model.getBoundary(model_dim_tags[0], oriented=False)
boundary_ids = [b[1] for b in boundary]
model.addPhysicalGroup(2, boundary_ids, tag=1)
# Add physical tag for the volume (Torus->tag=3, sphere->tag=2)
volume_entities = [model_dim_tags[0][0][1]]
model.addPhysicalGroup(3, volume_entities, tag=3)
volume_entities = [model_dim_tags[0][1][1]]
model.addPhysicalGroup(3, volume_entities, tag=2)
gmsh.option.setNumber('Mesh.MeshSizeMin', refinement)
gmsh.option.setNumber('Mesh.MeshSizeMax', refinement)
model.mesh.generate(dim=3)
### Solve the equation with the mesh ###
mesh, ct, ft = gmshio.model_to_mesh(model, MPI.COMM_WORLD, rank=0)
cdim = mesh.topology.dim
num_facets_owned_by_proc = mesh.topology.index_map(cdim).size_local
geometry_entitites = dolfinx.cpp.mesh.entities_to_geometry(mesh, cdim, np.arange(num_facets_owned_by_proc, dtype=np.int32), False)
points = mesh.geometry.x
Scalar_DG0 = FunctionSpace(mesh, ("DG", 0))
Vector_DG0 = VectorFunctionSpace(mesh, ("DG", 0))
dim3_tags = np.unique(ct.values)
mu = Function(Scalar_DG0)
J = Function(Vector_DG0)
# Define the current (J that should follow the spire) and the permeability µ (mu)
for tag in dim3_tags:
cells = ct.find(tag)
sum_mean_coord=[0,0,0]
if tag == 2 : #sphere vacuum
mu_ = 4 * np.pi*1e-7
J_ = 0
elif tag == 3: #torus
mu_ = 4 * np.pi*1e-7
J_ = 1
for num_cell in cells:
add_coord_x=[]
add_coord_y=[]
add_coord_z=[]
for node in geometry_entitites[num_cell]:
coord = points[node]
add_coord_x=coord[0]
add_coord_y=coord[1]
add_coord_z=coord[2]
mean_coord=[np.mean(add_coord_x),np.mean(add_coord_y),0]
sum_mean_coord = [sum(i) for i in zip(sum_mean_coord, mean_coord)]
r=np.linalg.norm(mean_coord)
rotate_coord=np.array([[0,-1,0],[1,0,0],[0,0,1]])@np.array(mean_coord)
J.x.array[num_cell*3:num_cell*3+3] = rotate_coord/r
mu.x.array[cells] = mu_
#Set the equation
V = FunctionSpace(mesh, (space, order))
boundary_facet = ft.find(1)
boundary_dofs = locate_dofs_topological(V, cdim-1, boundary_facet)
u0 = Function(V)
u0.x.set(0)
bc = dirichletbc(u0, boundary_dofs)
u = TrialFunction(V)
v = TestFunction(V)
a = ((1 / mu) * inner(rot(u), rot(v)) + 1e-6 * inner(u,v))* dx
a += 1e-6 * inner(u,v)*dx #had to add this term for N2curl, change nothing for N1curl
L = inner(J, v) * dx
A = Function(V)
problem = LinearProblem(a, L, u=A, bcs=[bc])
Ah = problem.solve()
### export in VTK ###
Int = VectorFunctionSpace(mesh, ("DG", order_interpolation))
A_fin = Function(Int)
A_fin.interpolate(Ah)
with io.VTKFile(proc, f"VTK/Ned_{mesh_name}_{space}{order}_DG{order_interpolation}A.pvd", "w") as VTK:
VTK.write_mesh(mesh)
VTK.write_function(A_fin)
B_exp = Expression(rot(Ah), Int.element.interpolation_points())
B = Function(Int)
B.interpolate(B_exp)
with io.VTKFile(proc, f"VTK/Ned_{mesh_name}_{space}{order}_DG{order_interpolation}B.pvd", "w") as VTK:
VTK.write_mesh(mesh)
VTK.write_function(B)