I meet a problem with implementing the unsteady NS equation with Newton method. I need to implement a periodical boundary condition on the inlet. However it seems that with Newton method, the convergence is not very well, the message shows that the SNES Function norm is always at a high value like XXe-1 to XXe+1 and it doesn’t go further down. Do you have any suggestions to improve the situation? Thank you!

I here enclose my code for checking, thank you! NewtonNS.edp (5.2 KB)

I cannot run your script because it needs an input file that you didn’t share. You mention that you have periodic boundary conditions, but I don’t see any periodic fespace in your file.

This is a txt file where I read in the inlet max velocity (it’s a parabolic profile). The way that I implemented the periodical BC is using this data repeatedly by + on(4, ux = uaugm[(tstep-1)%1000]*(y-1)^2-uaugm[(tstep-1)%1000], uy = 0);

Also due to the implementation of Newtons method, I calculate the augmentation of the velocity of each timestep from the velocity list ulist which is the argument uaugm here.

But I am not sure it is the correct implementation for the situation

Sorry I just realise that my question may not be precise. I may have used a wrong term. The periodic BC I mean is having a repeating pattern inlet velocity using the data like a pulsatile flow instead of having the outlet as another inlet.

Even the first nonlinear system cannot be solved. How did you implement your two varf? Does it work with a Newton linearization implemented “by hand”? My best advice would be to start from something that works, e.g., a steady-state Navier–Stokes script, and then add time discretization and then your forcing term. It’s difficult otherwise to debug these multiline varf.

I highly doubt your vRes varf is correct. Your Dirichlet boundary condition on label 4 does not depend on the previous residual, which can never go to 0, hence the stagnation. At convergence of SNESSolve, keep in mind that vRes(0, v) should evaluate to a vector whose norm is close to 0. Here, that’s impossible since you have on(4, ux = uaugm[(tstep-1)%1000]*(y-1)^2-uaugm[(tstep-1)%1000], uy = 0). Maybe it should be on(4, ux = ucx-uaugm[(tstep-1)%1000]*(y-1)^2-uaugm[(tstep-1)%1000], uy = 0)?