11 |
os.rmdir(dir_name) |
os.rmdir(dir_name) |
12 |
print "Restart files %s have been removed."%dir_name |
print "Restart files %s have been removed."%dir_name |
13 |
|
|
14 |
NE=20 |
|
15 |
|
import sys |
16 |
|
import time |
17 |
|
t1 = time.time() |
18 |
|
|
19 |
|
if (len(sys.argv)>=3): |
20 |
|
NE=int(sys.argv[2]) |
21 |
|
else: |
22 |
|
NE=20 |
23 |
|
|
24 |
|
if (len(sys.argv)>=2): |
25 |
|
solver=sys.argv[1] |
26 |
|
else: |
27 |
|
solver='PCG' |
28 |
|
|
29 |
|
if solver!='PCG': |
30 |
|
extratol=0.001 |
31 |
|
else: |
32 |
|
extratol=1 |
33 |
|
|
34 |
DIM=2 |
DIM=2 |
35 |
H=1. |
H=1. |
36 |
L=4*H |
L=4*H |
69 |
else: |
else: |
70 |
f=options.restart_dir |
f=options.restart_dir |
71 |
print "restart from directory ",f |
print "restart from directory ",f |
72 |
dom=LoadMesh("mesh.nc") |
try: |
73 |
|
dom=LoadMesh("mesh.nc") |
74 |
|
except: |
75 |
|
pass |
76 |
ff=open(os.path.join(f,"stamp.%d"%dom.getMPIRank()),"r").read().split(";") |
ff=open(os.path.join(f,"stamp.%d"%dom.getMPIRank()),"r").read().split(";") |
77 |
t=float(ff[0]) |
t=float(ff[0]) |
78 |
t_out=float(ff[1]) |
t_out=float(ff[1]) |
94 |
dom=Rectangle(int(ceil(L*NE/H)),NE,l0=L,l1=H,order=2, useFullElementOrder=True,optimize=True) |
dom=Rectangle(int(ceil(L*NE/H)),NE,l0=L,l1=H,order=2, useFullElementOrder=True,optimize=True) |
95 |
else: |
else: |
96 |
dom=Brick(int(ceil(L*NE/H)),int(ceil(L*NE/H)),NE,l0=L,l1=L,l2=H,order=2, useFullElementOrder=True,optimize=True) |
dom=Brick(int(ceil(L*NE/H)),int(ceil(L*NE/H)),NE,l0=L,l1=L,l2=H,order=2, useFullElementOrder=True,optimize=True) |
97 |
dom.dump("mesh.nc") |
try: |
98 |
|
dom.dump("mesh.nc") |
99 |
|
except: |
100 |
|
pass |
101 |
x=dom.getX() |
x=dom.getX() |
102 |
T=Scalar(1,ReducedSolution(dom)) |
T=Scalar(1,ReducedSolution(dom)) |
103 |
for d in range(DIM): |
for d in range(DIM): |
133 |
# set up velovity problem |
# set up velovity problem |
134 |
# |
# |
135 |
sp=StokesProblemCartesian(dom) |
sp=StokesProblemCartesian(dom) |
136 |
sp.setTolerance(TOL) |
sp.setTolerance(TOL*extratol) |
137 |
sp.setToleranceReductionFactor(TOL) |
sp.setToleranceReductionFactor(TOL) |
138 |
x2=ReducedSolution(dom).getX() |
x2=ReducedSolution(dom).getX() |
139 |
p=-RA*(x2[DIM-1]-0.5*x2[DIM-1]**2) |
p=-RA*(x2[DIM-1]-0.5*x2[DIM-1]**2) |
155 |
viscosity=exp(A*(1./(1+T.interpolate(Function(dom)))-2./3.)) |
viscosity=exp(A*(1./(1+T.interpolate(Function(dom)))-2./3.)) |
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') |
159 |
# v,p=sp.solve(v,p,show_details=VERBOSE, verbose=True,max_iter=500,solver='GMRES') |
#v,p=sp.solve(v,p,show_details=VERBOSE, verbose=True,max_iter=500,solver='GMRES') |
160 |
|
#v,p=sp.solve(v,p,show_details=VERBOSE, verbose=True,max_iter=500,solver='MINRES') |
161 |
|
v,p=sp.solve(v,p,show_details=VERBOSE, verbose=True,max_iter=500,solver=solver) |
162 |
|
|
163 |
for d in range(DIM): |
for d in range(DIM): |
164 |
print "range %d-velocity"%d,inf(v[d]),sup(v[d]) |
print "range %d-velocity"%d,inf(v[d]),sup(v[d]) |
210 |
else: |
else: |
211 |
file(os.path.join(new_restart_dir,"stamp.%d"%dom.getMPIRank()),"w").write("%e; %e; %s; %s; %e;\n"%(t, t_out, n, n_out, dt)) |
file(os.path.join(new_restart_dir,"stamp.%d"%dom.getMPIRank()),"w").write("%e; %e; %s; %s; %e;\n"%(t, t_out, n, n_out, dt)) |
212 |
removeRestartDirectory(old_restart_dir) |
removeRestartDirectory(old_restart_dir) |
213 |
|
elapsed = time.time() - t1 |
214 |
|
print "plot","\t",NE,"\t",elapsed |
215 |
|
|