I’m looking to solve navier-Stokes equation with an IMEX method which gives a convective term on the form

as this gives unconditional stability. I have tried to do this with a bilinear form which includes a gridfunction un. Below is is the setup of the spaces and convective part of the bilinear form:

V = VectorH1(mesh,order=k, dirichlet=“wall|cyl|inlet”)

Q = H1(mesh,order=k-1)

X = V*Q
(u,p), (v,q) = X.TnT()
gfu = GridFunction(X)
velocity, pressure = gfu.components
uin = CoefficientFunction((1.5*4

*y*(0.41-y)/(0.41

*0.41),0))*

velocity.Set(uin, definedon=mesh.Boundaries(“inlet”))

un = gfu.components[0]

conv = BilinearForm(X)

c = InnerProduct(Grad(u) * un, v)dx

velocity.Set(uin, definedon=mesh.Boundaries(“inlet”))

un = gfu.components[0]

conv = BilinearForm(X)

c = InnerProduct(Grad(u) * un, v)

conv += c

Assembly of conv and updating of un is done while time stepping. This method seems to work for quite a few time steps, until much of the data suddenly become nan values. Any help with figuring out how to do implement this scheme properly would be much appreciated!