Hello everybody,
I have a program where I have a boundary condition that depends on the value of one vector field. When I run as a single MPI process the program running (mpirun -np 1 python3 program.py) however if I run with mpirun -np N python3 program.py where N>1 the program not running. Does anybody have any idea why this happening?
Cod:
from dolfin import *
from fenics import *
from mshr import * #a racshoz kell
import numpy as np
import time
sadas=0
time_stp_counter=0
meroszam=[]
ido_mon=[]
v_max=[]
idolepes_cnt=0#Ez a plottolashoz kell
#uj parameterek
ct_par_t=5
T_futes=300+sadas*100
szam=0
szam1=1
########################## parameterek
alfa=1e-3
T=1#szimulacio ideje
ro_stp=1e23
N_stp=int(Tro_stp)
dt=1/ro_stp
l=20#szimulal terresz hoszusaga
d=20#szimulalt terresz magassaga
g=-9.81#gravitacios gyorsulas
T0=300
D=1e-23#Diffuzios allando
rho0=1
mu=1e-25
##########################
domain=Rectangle(Point(0,d/2),Point(l,-d/2))
i_max=10
mesh=generate_mesh(domain,128)#*2)
V=VectorFunctionSpace(mesh,‘P’,2)
Q=FunctionSpace(mesh,‘P’,1)
Qt=FunctionSpace(mesh,‘P’,2)
def boundary(x,on_boundary):
return on_boundary
pipe_and1=‘near(x[0],0)’#bemenet
pipe_and2=‘near(x[0],’+str(l)+’)’#kimenet
pipe_walls=‘near(x[1],’+str(-d/2)+’)||near(x[1],’+str(d/2)+’)’#falak sebesseg 0
bc_t_a1=‘300+’+str(T_futes)+’(-1/(pow(3,’+str(ct_par_t)+’(-x[1]+2))+1)+1/(pow(3,’+str(ct_par_t)+’(-x[1]-2))+1))’
#(-1/(pow(3,ct_par_t(-x+2))+1)+1/(pow(3,ct_par_t*(-x-2))+1))(-x+2)(x+2)0.1 // python szintaxis
bc_u_a1x=’(-1/(pow(3,’+str(ct_par_t)+’(-x[1]+2))+1)+1/(pow(3,’+str(ct_par_t)+’(-x[1]-2))+1))(-x[1]+2)(x[1]+2)0.4’+str(2.5)
bc_u_a1y=’(-1/(pow(3,’+str(ct_par_t)+’(-x[1]+2))+1)+1/(pow(3,’+str(ct_par_t)+’(-x[1]-2))+1))x[1]0.4’+str(4.5)
#-(x-d/2)(x+d/2)/1000.306916292607 //python szintaxis
bc_u_a2=’-(x[1]-’+str(d/2)+’)(x[1]+’+str(d/2)+’)/1000.306916292607’
#Hatarfeltetelek a sebessegre
bc_u2=DirichletBC(V,Expression((bc_u_a2,‘0’),degree=2),pipe_and2)
bc_u1=DirichletBC(V,Expression((bc_u_a1x,bc_u_a1y),degree=2),pipe_and1)
bc_uw=DirichletBC(V,Constant((0,0)),pipe_walls)
#Hatarfeltetelek a nyomasra
bc_p1=DirichletBC(Q,Constant(0),pipe_and1)
bc_p2=DirichletBC(Q,Constant(grho0l),pipe_and2)
#Homerseklet hatarfeltetel
bc_ta1=DirichletBC(Qt,Expression(bc_t_a1,degree=2),pipe_and1)
bc_tw=DirichletBC(Qt,Constant(T0),pipe_walls)
u=TrialFunction(V)
p=TrialFunction(Q)
T=TrialFunction(Qt)
v=TestFunction(V)
q=TestFunction(Q)
Tv=TestFunction(Qt)
u_n=Function(V)
u_s_top=Function(V)
p_n=Function(Q)
u_=Function(V)
p_=Function(Q)
T_n=Function(Qt)
Hoforras=Function(Qt)
f=Function(V)
n=FacetNormal(mesh)
k=Constant(dt)
mu=Constant(mu)
rho0=Constant(rho0)
#T0=Constant(T0)
D=Constant(D)
#Neuman_t=Constant(1e-3)#Homerseklet neuman hatarfeltetel
U=0.5*(u_n+u)
def epsilon(u):
return sym(nabla_grad(u))
def sigma(u,p):
return 2muepsilon(u)-p*Identity(len(u))
stt=‘0’#str(T_futes)+‘0.3(-1/(pow(3,’+str(ct_par_t)+’(-x[1]+2))+1)+1/(pow(3,’+str(ct_par_t)+’(-x[1]-2))+1))pow(3,-0.5x[0])’#‘0*pow(2.73,-pow(x[1]+0.3,2))’
stt=Expression(stt,degree=2)
Hoforras.assign(interpolate(stt,Qt))
def special_boundary_0(x,on_boundary):
return on_boundary and near(x[0],l) and u_n(x[0],x[1])[0]<0
def special_boundary_1(x,on_boundary):
return on_boundary and near(x[0],l) and ((u_n(x[0],x[1])[0]**2+u_n(x[0],x[1])[1]**2)>900)
bc_ts0=DirichletBC(Qt,Constant(T0),special_boundary_0)
bc_us0=DirichletBC(V,u_s_top,special_boundary_1)
bcu=[bc_u1,bc_uw,bc_us0]# Sebesseg hatarfeltetelek
bcp=[bc_p2]#nyomas hatarfeltetelek
bct=[bc_ta1,bc_tw,bc_ts0]#Homerseklet hatarfeltetelek
F1=rho0/(1+(T_n-T0)*alfa)*dot((u-u_n)/k,v)*dx
+rho0/(1+(T_n-T0)*alfa)*dot(dot(u_n,nabla_grad(u_n)),v)*dx
+inner(sigma(U,p_n),epsilon(v))*dx
-dot(rho0/(1+(T_n-T0)*alfa)*Constant((g,0)),v)dx
+dot(p_nn,v)ds-dot(munabla_grad(U)*n,v)*ds
a1=lhs(F1)
L1=rhs(F1)
a2=dot(nabla_grad§,nabla_grad(q))*dx
L2=dot(nabla_grad(p_n),nabla_grad(q))*dx-(1/k)*div(u_)qdx
a3=dot(u,v)*dx
L3=dot(u_,v)dx-kdot(nabla_grad(p_-p_n),v)*dx
F4=(((T-T_n)/k+dot(u_n,nabla_grad(T)))Tv+Ddot(nabla_grad(T),nabla_grad(Tv)))dx-HoforrasTvdx#-neuman_tTvds#-Ddot(nabla_grad(T),n)Tvds
a4=lhs(F4)
L4=rhs(F4)
#######a matrixok kiszamitasa elott uj erteket adok a homersekletnek
A1=assemble(a1)
A2=assemble(a2)
A3=assemble(a3)
A4=assemble(a4)
[bc.apply(A1) for bc in bcu]
[bc.apply(A2) for bc in bcp]
[bc.apply(A4) for bc in bct]
beta=5
terresz_fele=3
T=Function(Qt)
sc_t=‘300+’+str(T_futes)+’(-1/(pow(3,’+str(ct_par_t)+’(-x[1]+2))+1)+1/(pow(3,’+str(ct_par_t)+’*(-x[1]-2))+1))pow(3,-0.4x[0])’
sc_u=Constant((0,0))
#Homerseklet kezdeti feltetel
print(stt)
sc_t=Expression(sc_t,degree=2)
T_n.assign(interpolate(sc_t,Qt))
#Sebesseg kezdeti feltetel
u_n.assign(interpolate(sc_u,V))
u_.assign(interpolate(sc_u,V))
counter=0
‘’’
timeseries_u = TimeSeries(‘NS’+str(sadas)+’/velocity_series’)
timeseries_p = TimeSeries(‘NS’+str(sadas)+’/pressure_series’)
timeseries_T = TimeSeries(‘NS’+str(sadas)+’/Temperature_series’)
File(‘NS’+str(sadas)+’/mesh.xml.gz’) << mesh
‘’’
t=0
for n in range(10):
print(n)
t=t+dt
szam=szam+1
A4=assemble(a4)
[bc.apply(A4) for bc in bct]
b4 = assemble(L4)
[bc.apply(b4) for bc in bct]
solve(A4, T.vector(), b4, 'bicgstab', 'hypre_amg')
#solve(a4==L4,T)#,new_boundary_t)
T_n.assign(T)
A1=assemble(a1)
[bc.apply(A1) for bc in bcu]
b1 = assemble(L1)
[bc.apply(b1) for bc in bcu]
solve(A1, u_.vector(), b1, 'bicgstab', 'hypre_amg')
# Step 2: Pressure correction step
A2=assemble(a2)
[bc.apply(A2) for bc in bcp]
b2 = assemble(L2)
[bc.apply(b2) for bc in bcp]
solve(A2, p_.vector(), b2, 'bicgstab', 'hypre_amg')
# Step 3: Velocity correction step
A3=assemble(a3)
b3 = assemble(L3)
solve(A3, u_.vector(), b3, 'cg', 'sor')
adaptiv_meroszam=(u_-u_n)
adaptiv_meroszam=project(adaptiv_meroszam,V)
adaptiv_meroszam=adaptiv_meroszam.vector().max()
#Megj itt az is lehet a baj, hogy nem eleg suru a racs
#Kiprobalom ketszer olyan suru raccsal
meroszam.append(adaptiv_meroszam/dt)
ido_mon.append(dt)
v_max.append(u_n.vector().max())
u_n.assign(u_)
p_n.assign(p_)
if (u_n.vector().max()<30):
u_s_top.assign(u_n)
plot(T_n)