31 |
else: |
else: |
32 |
extratol=1 |
extratol=1 |
33 |
|
|
34 |
DIM=2 |
DIM=3 |
35 |
H=1. |
H=1. |
36 |
L=4*H |
L=4*H |
37 |
THETA=0.5 |
THETA=0.5 |
152 |
while t<T_END: |
while t<T_END: |
153 |
v_last=v*1 |
v_last=v*1 |
154 |
print "============== solve for v ========================" |
print "============== solve for v ========================" |
155 |
viscosity=exp(A*(1./(1+T.interpolate(Function(dom)))-2./3.)) |
viscosity=exp(A*(1./(1+T.interpolate(Function(dom)))-1./2.)) |
156 |
print "viscosity range :", inf(viscosity), sup(viscosity) |
print "viscosity range :", inf(viscosity), sup(viscosity) |
157 |
sp.initialize(f=T*(RA*unitVector(DIM-1,DIM)),eta=viscosity,fixed_u_mask=fixed_v_mask) |
sp.initialize(f=T*(RA*unitVector(DIM-1,DIM)),eta=viscosity,fixed_u_mask=fixed_v_mask) |
158 |
#v,p=sp.solve(v,p,show_details=VERBOSE, verbose=True,max_iter=500,solver='PCG') |
#v,p=sp.solve(v,p,show_details=VERBOSE, verbose=True,max_iter=500,solver='PCG') |
164 |
print "range %d-velocity"%d,inf(v[d]),sup(v[d]) |
print "range %d-velocity"%d,inf(v[d]),sup(v[d]) |
165 |
|
|
166 |
if t>=t_out: |
if t>=t_out: |
167 |
saveVTK("state.%d.vtu"%n_out,T=T,v=v,p=p) |
saveVTK("state.%d.vtu"%n_out,T=T,v=v) |
168 |
print "visualization file %d for time step %e generated."%(n_out,t) |
print "visualization file %d for time step %e generated."%(n_out,t) |
169 |
n_out+=1 |
n_out+=1 |
170 |
t_out+=DT_OUT |
t_out+=DT_OUT |