Hello,
Is there any tutorial where I can get my hands on that solves time dependent Navier Stokes equations using Taylor Hood pair (not with stepping), with DG discretization (not HDG) and where the nonlinear term is treated explicitly on the right hand side of the equation. I am trying to put the pieces together from the tutorials provided but there seems to be always a problem with the test and trial space definitions and the jump/average calculations. More specifically I have derived below-copied code from the tutorials:
>
> from ngsolve import *
>
> # viscosity
> nu = 0.001
>
> # timestepping parameters
> tau = 0.001
> tend = 10
>
> from netgen.geom2d import SplineGeometry
> geo = SplineGeometry()
> geo.AddRectangle( (0, 0), (2, 0.41), bcs = ("wall", "outlet", "wall", "inlet"))
> geo.AddCircle ( (0.2, 0.2), r=0.05, leftdomain=0, rightdomain=1, bc="cyl", maxh=0.02)
> mesh = Mesh( geo.GenerateMesh(maxh=0.07))
>
> mesh.Curve(3)
>
> V = VectorH1(mesh, order=2, dgjumps=True, dirichlet="wall|cyl|inlet")
> Q = H1(mesh,order=1)
>
> X=V*Q
>
> # Define test and trial functions for velocity and pressure
> u, v = V.TnT()
> p, q = Q.TnT()
>
> u = GridFunction(V)
> p = GridFunction(Q)
>
> jump_u = u-u.Other()
> jump_v = v-v.Other()
> n = specialcf.normal(2)
> mean_dudn = 0.5*n * (grad(u)+grad(u.Other()))
> mean_dvdn = 0.5*n * (grad(v)+grad(v.Other()))
>
> stokes = nu*InnerProduct(grad(u), grad(v))+div(u)*q+div(v)*p - 1e-10*p*q
> a = BilinearForm(X, symmetric=True)
> a += stokes*dx
> a.Assemble()
>
> # nothing here ...
> f = LinearForm(X)
> f.Assemble()
>
> # gridfunction for the solution
> gfu = GridFunction(X)
>
> # parabolic inflow at inlet:
> uin = CoefficientFunction( (1.5*4*y*(0.41-y)/(0.41*0.41), 0) )
> gfu.components[0].Set(uin, definedon=mesh.Boundaries("inlet"))
>
> # solve Stokes problem for initial conditions:
> inv_stokes = a.mat.Inverse(X.FreeDofs())
>
> res = f.vec.CreateVector()
> res.data = f.vec - a.mat*gfu.vec
> gfu.vec.data += inv_stokes * res
>
>
>
> # matrix for implicit Euler
> mstar = BilinearForm(X, symmetric=True)
>
> alpha = 4
> h = specialcf.mesh_size
> diffusion = +alpha*order**2/h*jump_u*jump_v * dx(skeleton=True) \
> +(-mean_dudn*jump_v-mean_dvdn*jump_u) * dx(skeleton=True) \
> +alpha*order**2/h*u*v * ds(skeleton=True) \
> +(-n*grad(u)*v-n*grad(v)*u)* ds(skeleton=True)
>
> mstar += tau*stokes*dx
> mstar += BilinearForm(diffusion).Assemble()
> mstar.Assemble()
> inv = mstar.mat.Inverse(X.FreeDofs(), inverse="sparsecholesky")
>
>
> # the non-linear term
> conv = BilinearForm(X, nonassemble = True)
> conv += (grad(u) * u) * v * dx
>
> # for visualization
> Draw (Norm(gfu.components[0]), mesh, "velocity", sd=3)
>
> # implicit Euler/explicit Euler splitting method:
> t = 0
> with TaskManager():
> while t < tend:
> print ("t=", t, end="\r")
>
> conv.Apply (gfu.vec, res)
> res.data += a.mat*gfu.vec
> gfu.vec.data -= tau * inv * res
>
> t = t + tau
> Redraw()
but I get the error:
Generate Mesh from spline geometry
Curve elements, order = 3
Traceback (most recent call last):
File "<string>", line 33, in <module>
File "/usr/lib/python3.10/dist-packages/ngsolve/utils.py", line 42, in grad
add = func.Operator("grad")
netgen.libngpy._meshing.NgException: Operator grad not overloaded for CF N5ngfem24OtherCoefficientFunctionE
Thank you for using NGSolve
Thank you in advance.