Hi,
I have trouble implementing a SIBC for an eddy current problem. Especially my code fails when I use the CGSolver. The error that I get looks like this:
ZGETRI::info = 1
ZGETRF::info = 1
I have made three options in my code: use SIBC, use reduced vector potential (A = A_0 + A_eddy), and use CGSolver. The options use SIBC does not work with use CGSolver. I hope that you can help to find my mistake.
The code:
from ngsolve import *
from netgen.occ import *
import numpy as np
import netgen.gui
use_SIBC = True
use_reduced_potential = True
use_CGSolver = True #False # using CGSolver is not working with SIBC. Why?
# inner sphere
sphere_in = Sphere(Pnt(0,0,0), 1)
sphere_in.mat("sphere inner")
sphere_in.solids.name = "sphere inner"
sphere_in.faces.name = "sphere inner"
sphere_in.faces.maxh = 0.3
# outer domain
air = Sphere(Pnt(0,0,0), 5)
air.solids.name = "air"
air.mat("air")
if use_SIBC:
geo = OCCGeometry((air - sphere_in)) # Subtract sphere from air domain (only solve in air)
else:
air = air - sphere_in
geo = OCCGeometry(Glue([air, sphere_in]))
# 2. Create Mesh and Define Finite Element Space
mesh = geo.GenerateMesh(maxh=1)
if use_SIBC == False:
mesh.BoundaryLayer(".*", [0.002, 0.004], material="sphere inner", domains="sphere inner") # add mesh layers
mesh.Curve(5)
mesh = Mesh(mesh)
#Draw(mesh, clipping = { "pnt" : (0,1,0), "vec" : (0,0,-1)})
#input("Press Enter to continue...")
# Use H(curl) elements for vector potential A
import math
fes = HCurl(mesh, order=2, nograds = False, complex=True)
A = fes.TrialFunction()
v = fes.TestFunction()
# Define Material Properties
mu_0 = 4 * np.pi * 1e-7 # Permeability of free space
sigma_in = 5.8e7 # Conductivity of copper
omega = 2 * np.pi * 50 # Frequency (50 Hz)
mu_r = { "air" : 1, "sphere inner": 1}
sigma = {"air": 1e-4, "sphere inner": sigma_in}
mu_0 = 4*math.pi*1e-7
nu_coef = [ 1/(mu_0*mu_r[mat]) for mat in mesh.GetMaterials() ]
sigma_coef = [sigma[mat] for mat in mesh.GetMaterials() ]
eta_coef = [ 1j*omega*sigma[mat] for mat in mesh.GetMaterials() ]
sigma = CoefficientFunction(sigma_coef)
nu = CoefficientFunction(nu_coef)
eta = CoefficientFunction(eta_coef)
# Define Bilinear Form (Weak Formulation) and the Right-Hand Side (Excitation Term)
A_s = CoefficientFunction((-y/2, x/2, 0)) # Impressed vector potential --> uniform B-field in z-direction: B_s=(0,0,1)
B_s = CoefficientFunction((0, 0, 1))
a = BilinearForm(fes, symmetric=False) #a = BilinearForm(fes, symmetric=True, condense=True)
f = LinearForm(fes)
if use_reduced_potential: # solving for reduced vector potential
if use_SIBC:
a += (1/mu_0) * ( curl(A) * curl(v) * dx + 1e-6*curl(A)*curl(v)*dx ) # Curl-curl term in air
a += sqrt(1j*omega*sigma_in/mu_0) * A.Trace() * v.Trace() * ds(definedon=mesh.Boundaries("sphere inner")) # SIBC on sphere surface
f += sqrt(1j*omega*sigma_in/mu_0) * A_s * v.Trace() * ds(definedon=mesh.Boundaries("sphere inner")) # Excitation term
else:
a += SymbolicBFI(eta*A*v + nu*curl(A)*curl(v) + 1e-6*nu*curl(A)*curl(v)) #a += eta*u*v*dx + nu*curl(u)*curl(v)*dx
f += eta*A_s*v*dx
else: # solving for full vector potential
if use_SIBC:
a += (1/mu_0) * ( curl(A) * curl(v) * dx + 1e-6*curl(A)*curl(v)*dx ) # Curl-curl term in air
a += sqrt(1j*omega*sigma_in/mu_0) * A.Trace() * v.Trace() * ds(definedon=mesh.Boundaries("sphere inner")) # SIBC on sphere surface
H_s = mesh.MaterialCF({}, default = B_s/mu_0) # External field
f += H_s*curl(v) * dx
else:
a += SymbolicBFI(eta*A*v + nu*curl(A)*curl(v) + 1e-6*nu*curl(A)*curl(v)) #a += eta*u*v*dx + nu*curl(u)*curl(v)*dx
H_s = mesh.MaterialCF({}, default = B_s/mu_0) # External field
f += H_s*curl(v) * dx
if use_CGSolver:
c = Preconditioner(a, type="bddc")
""" Preconditioners:
Jacobi (local) and block Jacobi
Direct solvers using sparse factorization (direct)
Geometric multigrid with different block-smoothers (multigrid)
Algebraic multigrid preconditioner (h1amg)
p-version element-level BDDC (bddc) """
# Solve the System
gfu = GridFunction(fes)
print("Solving...")
with TaskManager():
a.Assemble()
f.Assemble()
if use_CGSolver:
solver = CGSolver(mat=a.mat, pre=c.mat)
gfu.vec.data = solver * f.vec
else:
gfu.vec.data = a.mat.Inverse(freedofs=fes.FreeDofs()) * f.vec
#Draw(gfu, mesh, "Magnetic Vector Potential")
#input("Press Enter to continue...")
if use_reduced_potential:
B_eddy = curl(gfu)
else:
B_eddy = curl(gfu) - B_s
Draw(B_eddy.Norm(), mesh, 'B-field eddy')
input("Press Enter to continue...")
