Hello everyone,
I’m new using NGSolve (1 month) and I’m trying to solve a Taylor-Green Vortex problem in NGSolve. I don’t know if the problem is in the variational form or in the definition of the problem. I guess that it is in the boundary condition imposed. Please if someone can help me to see what it is, I’ll really appreciate it.
I’m using the exact solution to compare and compute the L2 Error of the numerical solution. Also, I’m using the same NS implicit/Explicit scheme shown in the tutorials: Navier Stokes Equation — NGS-Py 6.2.2404-18-gd2e5801de documentation (ngsolve.org)
Thank you in advance community 
I submit the code developed:
from ngsolve import *
import numpy as np
from ngsolve.webgui import Draw
import ipywidgets as widgets
from netgen.occ import *
from ngsolve.meshes import MakeStructured2DMesh
def F(t): return exp(-2 * np.pi**2 * t)
def vel(t):
u1 = sin(np.pi * x) * cos(np.pi * y) * F(t)
u2 = -sin(np.pi * y) * cos(np.pi * x) * F(t)
u_cond = CoefficientFunction((u1,u2))
return u_cond
def pre(t):
term1 = sin(2 * np.pi * x)**2
term2 = sin(2 * np.pi * y)**2
pre_cond = CoefficientFunction(-0.5 * (term1 + term2) * F(t)**2)
return pre_cond
#Mesh generation
L=1
bbname = ["Pressure_Point"]
#Geometry and Mesh parameters
n = 50; coord = [(0.1,0.1)]
#mesh = MakeStructured2DMesh(quads=True, nx=n,ny=n, mapping= lambda x,y : (L*x,L*y))
mesh = MakeStructured2DMesh(quads=True, nx=n,ny=n, mapping= lambda x,y : (L*x,L*y), bbpts=coord, bbnames=bbname)
V = VectorH1(mesh,order=2, dirichlet="left|bottom|right|top")
Q = H1(mesh,order=1)
X = V*Q
u,p = X.TrialFunction()
v,q = X.TestFunction()
nu = 1 # viscosity
stokes = (nu*InnerProduct(grad(u), grad(v))+ \
div(u)*q+div(v)*p - 1e-10*p*q)*dx
a = BilinearForm(stokes).Assemble()
#nothing here ...
f = LinearForm(X).Assemble()
#gridfunction for the solution
gfu = GridFunction(X)
#Boundary conditions
gfu.components[0].Set(vel(t), BND) # Velocity temporal-BC
gfu.components[1].Set(0, BBND) # Pressure Point Condition
#Initial Conditions
gfu.components[0].Set(vel(t))
tau = 0.01 # timestep
mstar = BilinearForm(u*v*dx + tau*stokes).Assemble()
inv = mstar.mat.Inverse(X.FreeDofs(), inverse="sparsecholesky")
conv = BilinearForm(X, nonassemble=True)
conv += (Grad(u) * u) * v * dx
t = 0; i = 0
tend = 1
scene = Draw (gfu.components[0], mesh)
scene_p = Draw (gfu.components[1], mesh)
tw = widgets.Text(value='t = 0')
display(tw)
wn = GridFunction(X)
p_exac = GridFunction(X)
scene_ex_v = Draw (p_exac.components[0])
scene_ex_p = Draw (p_exac.components[1])
u_n = gfu.components[0]
p_n = gfu.components[1]
with TaskManager():
while t < tend:
t = round(t,4)
res = conv.Apply(gfu.vec) + a.mat*gfu.vec
gfu.vec.data -= tau * inv * res
u_n1 = gfu.components[0]
p_n1 = gfu.components[1]
p_exac.components[0].Set(vel(t))
p_exac.components[1].Set(pre(t))
u_ex = p_exac.components[0]
p_ex = p_exac.components[1]
Vel_L2_Error = sqrt(Integrate((u_ex-u_n1)**2, mesh))
Pre_L2_Error = sqrt(Integrate((p_ex-p_n1)**2, mesh))
Vel_Rel_Error = sqrt(Integrate((u_n1-u_n)**2, mesh))
Pre_Rel_Error = sqrt(Integrate((p_n1-p_n)**2, mesh))
print('Velocity L2 Error:',Vel_L2_Error)
print('Pressure L2 Error:',Pre_L2_Error)
wn.components[0].Set(gfu.components[0])
wn.components[1].Set(gfu.components[1])
u_n = wn.components[0]
p_n = wn.components[1]
t = t + tau; i = i + 1
if i%10 == 0:
scene.Redraw()
scene_p.Redraw()
scene_ex_p.Redraw()
scene_ex_v.Redraw()
if i%10 == 0: tw.value = "t = " + str(t)
tw.value = "t = "+str(tend)