Problem with .Assemble()

from ngsolve import * from netgen import NgOCC from math import * import sys geo = NgOCC.LoadOCCGeometry("Normale_Geometrie_3.step") mesh= Mesh(geo.GenerateMesh()) dirichletstring = "bc_0|bc_1|bc_2" for i in range(4,103): dirichletstring = dirichletstring + "|" + "bc_" + "%d" % (i) V = VectorH1(mesh,order=3, dirichlet=dirichletstring) Q = H1(mesh,order=2) Q_alt = H1(mesh, order=2, dirichlet = dirichletstring) N = NumberSpace(mesh) X = FESpace([V,Q]) #constants density = 0.006984526006623 sigma_k = 1 sigma_e = 1.3 C_e1 = 1.44 C_e2 = 1.92 C_mu = 0.09 #constant viscosity nu = 35.1524056250947 #relaxation parameters omega_k = 0.4 omega_T = 0.4 omega_e = 0.4 #trial and test functions u,p = X.TrialFunction() v,q = X.TestFunction() u_star = V.TrialFunction() phi = Q.TrialFunction() l= V.TrialFunction() k_star = Q_alt.TrialFunction() e_star = Q_alt.TrialFunction() e = Q_alt.TrialFunction() v_new = Q_alt.TestFunction() #gridfunctions u_star_gfu = GridFunction(V) phi_gfu = GridFunction(Q) l_gfu = GridFunction(V) k_star_gfu = GridFunction(Q_alt) k_gfu_last = GridFunction(Q_alt) k_gfu_next = GridFunction(Q_alt) k_iter_gfu_last = GridFunction(Q_alt) k_iter_gfu_next = GridFunction(Q_alt) e_star_gfu = GridFunction(Q_alt) e_gfu_last = GridFunction(Q_alt) e_gfu_next = GridFunction(Q_alt) e_iter_gfu_last = GridFunction(Q_alt) e_iter_gfu_next = GridFunction(Q_alt) p_k_gfu = GridFunction(Q_alt) f_2_gfu = GridFunction(Q_alt) f_mu_gfu = GridFunction(Q_alt) Re_T_gfu = GridFunction(Q_alt) gfu_last = GridFunction(X) gfu_next = GridFunction(X) nu_T_star_gfu = GridFunction(N) nu_T_gfu_last = GridFunction(N) nu_T_gfu_next = GridFunction(N) nu_T_iter_gfu_last = GridFunction(N) nu_T_iter_gfu_next = GridFunction(N) #boundary conditions u_inlet = CoefficientFunction( (0,0,5)) gfu_last.components[0].Set(u_inlet, definedon=mesh.Boundaries("bc_103")) k_inlet = 0.01 * 25 e_inlet = C_mu * sqrt(k_inlet) * k_inlet / 1100 k_star_gfu.Set(CoefficientFunction(k_inlet), definedon = mesh.Boundaries("bc_103")) k_gfu_last.Set(CoefficientFunction(k_inlet), definedon = mesh.Boundaries("bc_103")) k_gfu_next.Set(CoefficientFunction(k_inlet), definedon = mesh.Boundaries("bc_103")) k_iter_gfu_last.Set(CoefficientFunction(k_inlet), definedon = mesh.Boundaries("bc_103")) k_iter_gfu_next.Set(CoefficientFunction(k_inlet), definedon = mesh.Boundaries("bc_103")) e_star_gfu.Set(CoefficientFunction(e_inlet), definedon = mesh.Boundaries("bc_103")) e_gfu_last.Set(CoefficientFunction(e_inlet), definedon = mesh.Boundaries("bc_103")) e_gfu_next.Set(CoefficientFunction(e_inlet), definedon = mesh.Boundaries("bc_103")) e_iter_gfu_last.Set(CoefficientFunction(e_inlet), definedon = mesh.Boundaries("bc_103")) e_iter_gfu_next.Set(CoefficientFunction(e_inlet), definedon = mesh.Boundaries("bc_103")) #initial conditions mixing_length = 1100 k_0 = (nu/mixing_length)* (nu/mixing_length) e_0 = C_mu * sqrt(k_0) * k_0 / mixing_length Re_T_0 = k_0 * k_0 / (nu * e_0) f_mu_0 = exp(-3.4/((1+Re_T_0/50)*(1+Re_T_0/50))) k_gfu_last.Set(CoefficientFunction(k_0)) k_iter_gfu_last.Set(CoefficientFunction(k_0)) k_iter_gfu_next.Set(CoefficientFunction(k_0)) e_gfu_last.Set(CoefficientFunction(k_0)) e_iter_gfu_last.Set(CoefficientFunction(k_0)) e_iter_gfu_next.Set(CoefficientFunction(k_0)) nu_T_gfu_last.Set(CoefficientFunction(C_mu * f_mu_0 * (k_0 * k_0)/e_gfu_last)) #time specification time= 0.0 dt= 0.1 timeend = 0.5 Draw(gfu_last.components[0],mesh,"u") while time < timeend: print("time = ", time) #equation for u_star with beta=0 u_star_transient_i = InnerProduct(u_star,v) a_u_star = BilinearForm(V) a_u_star += (u_star_transient_i + InnerProduct(grad(u_star),grad(v)))*dx a_u_star.Assemble()