- Thank you received: 0
Problem with .Assemble()
3 years 2 months ago #3978
by leon123
Problem with .Assemble() was created by leon123
Hello everyone,
I am trying to implement a k-epsilon model for a turbulent flow in a three dimensional geometry.
When I try to run code, the program runs up to an assemble statement and then just exits without giving me an error message. Does someone have any insights on why the program just breaks doing an assemble statement.
Here is my code.
At this assemble statement the program just exits and I don't know why.
Thanks for your help.
Leon
I am trying to implement a k-epsilon model for a turbulent flow in a three dimensional geometry.
When I try to run code, the program runs up to an assemble statement and then just exits without giving me an error message. Does someone have any insights on why the program just breaks doing an assemble statement.
Here is my code.
Code:
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()
Thanks for your help.
Leon
Time to create page: 0.110 seconds