Hi there,
As the title suggests, I am working on a 1D problem in a 2D domain with interface conditions in cylindrical coordinates. I suspect there is a problem with the implementation of the interface conditions, and I would be very grateful if anyone could help me. I have looked through the ufl documentation and the discourse quite a lot, but I am still stuck.
This is the domain
If A_1 is the solution on \Omega_{left}, A_2 is the solution on \Omega_{middle}, and A_3 is the solution on \Omega_{right}, I is a positive real constant and C_1 is a complex constant, then the PDEs are
-\frac{1}{\mu_1}\frac{d}{dr}\left(\frac{1}{r}\frac{d(rA_1)}{dr}\right) + i \omega \sigma_1 A_1 = 0 \qquad \text{ if } r_0 < r< r_1 \\
-\frac{1}{\mu_0}\frac{d}{dr}\left(\frac{1}{r}\frac{d(rA_2)}{dr}\right) = 0\qquad \text{ if } r_1 < r < r_2 \\
-\frac{1}{\mu_0}\frac{d}{dr}\left(\frac{1}{r}\frac{d(rA_3)}{dr}\right) = 0\qquad \text{ if } r_2 < r < r_3 \\
A _1= 0 \qquad \text{ on } \Gamma_{left}\\
A_3 = C_1 \qquad \text{ on } \Gamma_{right} \\
A_1 = A_2 \qquad \text{ on } \Gamma_{left \ interior} \\
A_2 = A_3 \qquad \text{ on } \Gamma_{right\ interior} \\
\frac{1}{\mu_1}\frac{1}{r}\frac{d(rA_1)}{dr} = \frac{1}{\mu_0}\frac{1}{r}\frac{d(rA_2)}{dr} \qquad \text{ on } \Gamma_{left \ interior} \\
\frac{1}{\mu_0}\frac{1}{r}\left(\frac{d(rA_2)}{dr} -\frac{d(rA_3)}{dr}\right) = I \qquad \text{ on } \Gamma_{right\ interior}
A minimum working example of the code is here:
import dolfinx
import gmsh
from dolfinx import io, fem
from mpi4py import MPI
from petsc4py.PETSc import ScalarType
from ufl import Measure
import numpy as np
from helpers import functions
import meshio
from scipy import special as bessel
from ufl import inner, SpatialCoordinate, FacetNormal, TrialFunction, TestFunction
functions.clear_output()
gmsh.initialize()
lc = 0.00625
r1 = 0.1
r2 = 0.2
r3 = 0.8
length = 0.5
## Geometry
geo = gmsh.model.geo
geo.addPoint(0, length, 0, lc, 1)
geo.addPoint(r1, length, 0, lc, 2)
geo.addPoint(r2, length, 0, lc, 3)
geo.addPoint(r3, length, 0, lc, 4)
geo.addPoint(0, 0, 0, lc, 5)
geo.addPoint(r1, 0, 0, lc, 6)
geo.addPoint(r2, 0, 0, lc, 7)
geo.addPoint(r3, 0, 0, lc, 8)
geo.addLine(1, 2, 1) # top horizontal lines
geo.addLine(2, 3, 2)
geo.addLine(3, 4, 3)
geo.addLine(5, 6, 4) # bottom horizontal lines
geo.addLine(6, 7, 5)
geo.addLine(7, 8, 6)
geo.addLine(1, 5, 7) # vertical lines downward
geo.addLine(2, 6, 8)
geo.addLine(3, 7, 9)
geo.addLine(4, 8, 10)
geo.addCurveLoop([1, 8, -4, -7], 1)
geo.addCurveLoop([2, 9, -5, -8], 2)
geo.addCurveLoop([3, 10, -6, -9], 3)
geo.addPlaneSurface([1], 1)
geo.addPlaneSurface([2], 2)
geo.addPlaneSurface([3], 3)
gmsh.model.geo.synchronize()
left_boundary = 1
left_interface = 2
right_interface = 3
right_boundary = 4
bottom_boundary = 5
top_boundary = 6
left_area = 1
middle_and_right_area = 2
gmsh.model.addPhysicalGroup(1, [7], left_boundary)
gmsh.model.addPhysicalGroup(1, [8], left_interface)
gmsh.model.addPhysicalGroup(1, [9], right_interface)
gmsh.model.addPhysicalGroup(1, [10], right_boundary)
gmsh.model.addPhysicalGroup(1, [1,2,3], top_boundary)
gmsh.model.addPhysicalGroup(1, [4,5,6], bottom_boundary)
gmsh.model.addPhysicalGroup(2, [1], left_area)
gmsh.model.addPhysicalGroup(2, [2, 3], middle_and_right_area)
gmsh.model.mesh.generate(2)
gmsh.write(f"../output/mwe_mesh.msh")
gmsh.open(f"../output/mwe_mesh.msh")
gmsh.model.setCurrent("mwe_mesh")
print()
print("Model " + gmsh.model.getCurrent() + " (" + str(gmsh.model.getDimension()) + "D)")
print()
def create_mesh(mesh, cell_type, prune_z=False):
cells = mesh.get_cells_type(cell_type)
cell_data = mesh.get_cell_data("gmsh:physical", cell_type)
points = mesh.points[:, :2] if prune_z else mesh.points
out_mesh = meshio.Mesh(
points=points, cells={cell_type: cells}, cell_data={"name_to_read": [cell_data]}
)
return out_mesh
msh = meshio.read(f"../output/mwe_mesh.msh") # read into meshio
# Create and save one file for the mesh, and one file for the facets
triangle_mesh = create_mesh(msh, "triangle", prune_z=True) # prune = True because it is a 2D mesh
line_mesh = create_mesh(msh, "line", prune_z=True)
meshio.write("../output/mesh.xdmf", triangle_mesh) # write the mesh
meshio.write("../output/mt.xdmf", line_mesh)
with io.XDMFFile(MPI.COMM_WORLD, "../output/mesh.xdmf", "r") as xdmf:
domain = xdmf.read_mesh(name="Grid")
ct = xdmf.read_meshtags(domain, name="Grid") # these are the meshtags for the different 2D physical groups!
domain.topology.create_connectivity(domain.topology.dim, domain.topology.dim-1)
# writing so that we can debug in paraview
with io.XDMFFile(MPI.COMM_WORLD, "../output/mwe_cell_tags.xdmf", "w") as xdmf:
xdmf.write_mesh(domain)
xdmf.write_meshtags(ct, domain.geometry)
domain.topology.create_connectivity(domain.topology.dim-1 , domain.topology.dim)
with io.XDMFFile(MPI.COMM_WORLD, "../output/mt.xdmf", "r") as xdmf:
ft = xdmf.read_meshtags(domain, name="Grid") # these are the meshtags for the different 1D physical groups!
with io.XDMFFile(MPI.COMM_WORLD, "../output/mwe_facet_tags.xdmf", "w") as xdmf:
xdmf.write_mesh(domain)
xdmf.write_meshtags(ft, domain.geometry) # these are the meshtags for the different 1D physical groups!
V = fem.functionspace(domain, ("CG", 1))
interfaces = Measure("dS", domain=domain, subdomain_data=ft)
subdomains = Measure("dx", domain=domain, subdomain_data=ct)
left_facets = ft.find(left_boundary)
left_dofs = fem.locate_dofs_topological(V, domain.topology.dim - 1, left_facets)
right_facets = ft.find(right_boundary)
right_dofs = fem.locate_dofs_topological(V, domain.topology.dim - 1, right_facets)
mu0 = 4 * np.pi * 1e-7 # constants for the solution
mu1 = 259.47 * mu0
sigma1 = 5.911563017e6
omega = 37000
current = 30460
bcs = [fem.dirichletbc(ScalarType(0 + 0 * 1j), left_dofs, V), fem.dirichletbc(ScalarType(0.0007180974632771937-4.002859810663115e-07j), right_dofs, V)]
I = fem.Function(V)
I.x.array[:] = np.full_like(I.x.array, current)
x = SpatialCoordinate(domain)
n = FacetNormal(domain)
A = TrialFunction(V)
v = TestFunction(V)
dAdx = 1 / x[0] * (x[0] * A).dx(0)
dvdx = 1 / x[0] * (x[0] * v).dx(0)
# the weak form of the PDE with cylindrical coordinates
a = inner(1 / mu1 * dAdx, dvdx) * x[0] * subdomains(left_area) \
+ inner(1 / mu0 * dAdx, dvdx) * x[0] * subdomains(middle_and_right_area) \
+ 1j * omega * sigma1 * inner(A, v) * x[0] * subdomains(left_area) \
+ inner(inner(1 / mu1 * dAdx('+'), n[0]('+')), v('+')) * interfaces(left_interface) \
+ inner(inner(1 / mu0 * dAdx('-'), n[0]('-')), v('+')) * interfaces(left_interface) \
+ inner(inner(1 / mu0 * dAdx('+'), n[0]('+')), v('+')) * interfaces(right_interface)\
+ inner(inner(1 / mu0 * dAdx('-'), n[0]('-')), v('+')) * interfaces(right_interface) \
L = inner(I("+"), v("+")) * interfaces(right_interface)
problem_A = fem.petsc.LinearProblem(
a,
L,
bcs=bcs,
petsc_options={"ksp_type": "preonly", "pc_type": "lu"},
)
A_new = problem_A.solve()
xdmf = io.XDMFFile(domain.comm, "../output/mwe_A.xdmf", "w", encoding=io.XDMFFile.Encoding.ASCII)
xdmf.write_mesh(domain)
xdmf.write_function(A_new, 0)
xdmf.close()
The problem is that the solution this produces doesn’t match the analytical solution, and most importantly, does not converge to the analytical solution. See pic with analytic on the left, numerical on the right: the numerical does not reach the same maximum as the analytic.
I guess my question is: is this the correct way to implement the interface conditions? I think that might be where I am going wrong. Thank you for your time!
Best wishes,
Katie