Hello,
I am solving the PDE
\partial_t m - \nabla^2 m = -ifGm - m
where m=m_x + i m_y (i^2=-1 used just to simplify the notation), f=f(t) usually a pulse function and G(\vec{x})=\vec{g}\cdot\vec{x}. With initial condition and Robin boundary conditions
m(t=0) = 1
(\nabla m\cdot\hat{n} + m)|_{\vec{x}\in\partial \Omega} = 0
I based my time-stepping on the Cahn-Hilliard demo which subclasses NonlinearProblem and it works as it should as it allows to update the value of f(t) as the time-stepping goes on.
The issue I would appreciate some help with is that my code is not working for a type of mesh that I need while it does work for other meshes like UnitSquareMesh or even more complicated ones
The MWE:
import dolfin as fem
import numpy as np
import matplotlib.pyplot as plt
def BT(mesh):
# time-stepping params
dt = 0.001
tf = 0.200
nt = int(round(tf/dt)) + 1
t_arr = np.linspace(0., tf, nt)
theta = 0.5 # theta-method; 0.5 for Crank-Nicolson
# the function space needed for the problem
V = fem.FunctionSpace(mesh, fem.MixedElement(
(fem.FiniteElement('P', mesh.ufl_cell(), 1),
fem.FiniteElement('P', mesh.ufl_cell(), 1),
)
)
)
# define test functions
u, v = fem.TestFunctions(V)
# the solution at the next and current time-step:
m = fem.Function(V, name= 'm' )
mn = fem.Function(V, name= 'mn')
# set initial condition: m= (1., 0.) everywhere
m_local = m.vector().get_local()
m_local[::2] = 1.
m.vector().set_local(m_local)
m.vector().apply('insert')
# split the function into scalar-valued functions
mx, my = fem.split(m)
mxn, myn = fem.split(mn)
# for using the theta-method, define
mxnt = theta*mx + (1. - theta)*mxn
mynt = theta*my + (1. - theta)*myn
# we want a time- and position-dependent coefficient
fG = fem.Expression('f*G*(gx*x[0] + gy*x[1])',
f = 1., # for simplicity assume f(t) = 1. for all t
G = 1., # a low value
gx = 0.,
gy = 1.,
domain = mesh,
degree = 1,
)
# the variational form
L1 = (1./dt)*(mx-mxn)*u*fem.dx - fG*mynt*u*fem.dx + mxnt*u*fem.dx + \
fem.dot(fem.grad(mxnt), fem.grad(u))*fem.dx + mxnt*u*fem.ds
L2 = (1./dt)*(my-myn)*v*fem.dx + fG*mxnt*v*fem.dx + mynt*v*fem.dx + \
fem.dot(fem.grad(mynt), fem.grad(v))*fem.dx + mynt*v*fem.ds
L = L1 + L2
# problem setup, inspired by the Cahn-Hilliard demo
class Interface(fem.NonlinearProblem):
def __init__(self, a, L):
self.a = a
self.L = L
super().__init__()
def F(self, b, x):
fem.assemble(self.L, tensor= b)
def J(self, A, x):
fem.assemble(self.a, tensor= A)
dm = fem.TrialFunction(V)
a = fem.derivative(L, m, dm)
problem = Interface(a, L)
solver = fem.NewtonSolver()
solver.parameters['linear_solver'] = 'lu'
solver.parameters['convergence_criterion'] = 'incremental'
solver.parameters['relative_tolerance'] = 1.e-6
# time-stepping
mt = np.zeros(nt); mt[0] = fem.assemble(fem.sqrt(m[0]*m[0]+m[1]*m[1])*fem.dx)
for i, t in enumerate(t_arr[1:]):
mn.vector()[:] = m.vector()
solver.solve(problem, m.vector())
mt[i+1] = fem.assemble(fem.sqrt(m[0]*m[0]+m[1]*m[1])*fem.dx)
# plot results
fig, ax = plt.subplots()
plt.xlabel('t')
plt.ylabel('m(t)')
plt.plot(t_arr, mt/mt[0], label= 'calculated')
plt.legend()
bmesh = fem.BoundaryMesh(mesh, 'exterior')
area = 0.
for cell in fem.cells(bmesh):
area += cell.volume()
volume = 0.
for cell in fem.cells(mesh):
volume += cell.volume()
T2s = 1./(1. + (area/volume))
ax_error = ax.twinx()
ax_error.set_ylabel('error')
ax_error.plot(t_arr, np.abs((mt/mt[0])-np.exp(-t_arr/T2s)), 'r', label= 'error')
plt.legend(loc= (0.8, 0.8))
plt.show()
# create mesh
mesh = fem.UnitSquareMesh(20, 20)
# scale the mesh
mesh.scale(1.e-1)
BT(mesh) # this works
which results in
whereas, using the mesh I created with gmsh (download at https://www.dropbox.com/s/hmpknpafd1ygws7/mesh.zip?dl=0) I get:
# load the mesh
mesh = fem.Mesh()
with fem.XDMFFile('mesh.xdmf') as infile:
infile.read(mesh)
# scale the mesh
mesh.scale(1.e-2)
# plot for reference
fem.plot(mesh, linewidth= 0.1)
# solve
BT(mesh)
results in the plot of the mesh (for reference here)
and the error
RuntimeError:
*** -------------------------------------------------------------------------
*** DOLFIN encountered an error. If you are not able to resolve this issue
*** using the information listed below, you can ask for help at
***
*** fenics-support@googlegroups.com
***
*** Remember to include the error message listed below and, if possible,
*** include a *minimal* running example to reproduce the error.
***
*** -------------------------------------------------------------------------
*** Error: Unable to solve nonlinear system with NewtonSolver.
*** Reason: Newton solver did not converge because maximum number of iterations reached.
*** Where: This error was encountered inside NewtonSolver.cpp.
*** Process: 0
***
*** DOLFIN version: 2019.1.0
*** Git changeset: unknown
*** -------------------------------------------------------------------------
I can’t see anything obviously wrong with this mesh, I can even use it for simpler problems not involving time-stepping, i.e., non-homogenous Helmholtz equation:
\nabla^2u-u=1
with Robin B.C.
(\nabla u \cdot \hat{n} + u)|_{\vec{x}\in\partial\Omega} = 0
# load the mesh
mesh = fem.Mesh()
with fem.XDMFFile('/tmp/mesh.xdmf') as infile:
infile.read(mesh)
# scale the mesh
mesh.scale(1.e-2)
# scalar function space
V = fem.FunctionSpace(mesh, 'CG', 1)
# trial and test functions
u = fem.TrialFunction(V)
v = fem.TestFunction(V)
# variational form
a = (fem.dot(fem.grad(u), fem.grad(v)) + u*v)*fem.dx + u*v*fem.ds; L = v*fem.dx
# placeholder for solution
m = fem.Function(V)
# solve
fem.solve(a==L, m)
# view solution
fem.plot(m, levels= 256)
I would appreciate any help with this
Cheers