The solution diverges with the following nonlinear problem.
Are there any techniques to make it converge?
from numpy import *
mu0 = 4*pi*1e-7
aa = 0.002
sig = 1e7
Np = 101
time = linspace(0.00, 0.02, 4000)
dt = time[1]-time[0]
from netgen.occ import *
from netgen.webgui import Draw as DrawGeo
Wp = WorkPlane()
Wp.MoveTo(0, 0).Circle(aa)
geo = Wp.Face()
geo.name = 'Circle'
geo.edges.name = 'Circle'
geo.maxh = 0.1
geo.vertices.Nearest((1, 0)).name = "gnd"
import ngsolve
from ngsolve import *
from ngsolve.webgui import Draw
mesh = Mesh(OCCGeometry(geo, dim=2).GenerateMesh(meshsize.very_fine)).Curve(3)
print(mesh.GetBoundaries())
print(mesh.GetMaterials())
fes = H1(mesh, order=2, dirichlet_bbnd='gnd')
A = fes.TrialFunction()
W = fes.TestFunction()
gfA = GridFunction(fes)
gfA.Set(1)
Vs = lambda t: 0.5*tanh((5/aa) * (t - 0.004)) + 0.5
By = []
for t in time:
f_ = CoefficientFunction(-((2.0*Vs(t))/(500*mu0)+0.5*(2.0*Vs(t))**13))
a = BilinearForm(fes, symmetric=True)
a += SymbolicBFI(grad(W)*(1/(500*mu0)+0.5*(Norm(grad(A))**12))*grad(A))
a += SymbolicBFI(W*(2*sig/dt)*A)
a += - W * f_ * x/aa * ds(definedon=mesh.Boundaries('Circle'))
c = Preconditioner(a, type="bddc")
gfAh = GridFunction(fes)
gfAh.Set(gfA)
res = gfA.vec.CreateVector()
dA = gfA.vec.CreateVector()
for niter in range(100):
a.Apply(gfAh.vec, res)
a.AssembleLinearization(gfAh.vec)
solvers.CG(sol=dA.data, rhs=res.data, mat=a.mat, pre=c.mat, printrates='', tol=1e-10, maxsteps=10000)
# dA.data = a.mat.Inverse(fes.FreeDofs()) * res
gfAh.vec.data -= dA
if Norm(dA) < 1e-10:
print(f't = {t} niter = {niter} err = {Norm(dA):.3e}')
break
gfAh = 2*gfAh - gfA
gfA = GridFunction(fes)
gfA.Set(gfAh)
gfAh = GridFunction(fes)
gfAh.Set(gfA)
gfBy = -grad(gfA)[0]
x_ = linspace(0, aa, Np)
By.append([ gfBy(mesh(x,0)) for x in x_ if mesh.Contains(x)])