| 1 |
from esys.escript import * |
| 2 |
from esys.escript.models import TemperatureCartesian, StokesProblemCartesian |
| 3 |
from esys.finley import Rectangle, Brick, LoadMesh |
| 4 |
from optparse import OptionParser |
| 5 |
from math import pi, ceil |
| 6 |
def removeRestartDirectory(dir_name): |
| 7 |
if os.path.isdir(dir_name): |
| 8 |
for root, dirs, files in os.walk(dir_name, topdown=False): |
| 9 |
for name in files: os.remove(os.path.join(root,name)) |
| 10 |
for name in dirs: os.remove(os.path.join(root,name)) |
| 11 |
os.rmdir(dir_name) |
| 12 |
print "Restart files %s have been removed."%dir_name |
| 13 |
|
| 14 |
|
| 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=3 |
| 35 |
H=1. |
| 36 |
L=4*H |
| 37 |
THETA=0.5 |
| 38 |
TOL=1.e-3 |
| 39 |
PERTURBATION=0.1 |
| 40 |
T_END=0.3 |
| 41 |
DT_OUT=T_END/500 |
| 42 |
VERBOSE=False |
| 43 |
RA=1.e5 # Rayleigh number |
| 44 |
A=22. # Arenious number |
| 45 |
DI = 0. # dissipation number |
| 46 |
SUPG=False |
| 47 |
create_restartfiles_every_step=10 |
| 48 |
|
| 49 |
print "total number of elements = ",NE**DIM*int(L/H)**(DIM-1) |
| 50 |
|
| 51 |
# read options: |
| 52 |
parser = OptionParser(usage="%prog [-r [DIR]]") |
| 53 |
parser.add_option("-r", "--restart", dest="restart", help="restart from latest directory. It will be deleted after a new directory has been created.", default=False, action="store_true") |
| 54 |
parser.add_option("-d", "--dir", dest="restart_dir", help="restart from directory DIR. The directory will not be deleted but new restart directories are created.",metavar="DIR", default=None) |
| 55 |
(options, args) = parser.parse_args() |
| 56 |
restart=options.restart or (options.restart_dir !=None) |
| 57 |
# |
| 58 |
# set up domain: |
| 59 |
# |
| 60 |
if restart: |
| 61 |
if options.restart_dir ==None: |
| 62 |
restart_files=[] |
| 63 |
for f in os.listdir("."): |
| 64 |
if f.startswith("restart"): restart_files.append(f) |
| 65 |
if len(restart_files)==0: |
| 66 |
raise IOError,"no restart files" |
| 67 |
restart_files.sort() |
| 68 |
f=restart_files[-1] |
| 69 |
else: |
| 70 |
f=options.restart_dir |
| 71 |
print "restart from directory ",f |
| 72 |
try: |
| 73 |
dom=LoadMesh("mesh.nc") |
| 74 |
except: |
| 75 |
pass |
| 76 |
ff=open(os.path.join(f,"stamp.%d"%dom.getMPIRank()),"r").read().split(";") |
| 77 |
t=float(ff[0]) |
| 78 |
t_out=float(ff[1]) |
| 79 |
n=int(ff[2]) |
| 80 |
n_out=int(ff[3]) |
| 81 |
dt=float(ff[4]) |
| 82 |
v=load(os.path.join(f,"v.nc"),dom) |
| 83 |
p=load(os.path.join(f,"p.nc"),dom) |
| 84 |
T=load(os.path.join(f,"T.nc"),dom) |
| 85 |
if n>1: |
| 86 |
dt_a=float(ff[5]) |
| 87 |
a=load(os.path.join(f,"a.nc"),dom) |
| 88 |
else: |
| 89 |
dt_a=None |
| 90 |
a=None |
| 91 |
if dom.getMPIRank()==0: nusselt_file=open("nusselt.csv","a") |
| 92 |
else: |
| 93 |
if DIM==2: |
| 94 |
dom=Rectangle(int(ceil(L*NE/H)),NE,l0=L,l1=H,order=2, useFullElementOrder=True,optimize=True) |
| 95 |
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) |
| 97 |
try: |
| 98 |
dom.dump("mesh.nc") |
| 99 |
except: |
| 100 |
pass |
| 101 |
x=dom.getX() |
| 102 |
T=Scalar(1,ReducedSolution(dom)) |
| 103 |
for d in range(DIM): |
| 104 |
if d == DIM-1: |
| 105 |
T*=sin(x[d]/H*pi) |
| 106 |
else: |
| 107 |
T*=cos(x[d]*((d+1)/L*pi)) |
| 108 |
|
| 109 |
T=1.-x[DIM-1]+PERTURBATION*T |
| 110 |
v=Vector(0,Solution(dom)) |
| 111 |
if dom.getMPIRank() ==0: nusselt_file=open("nusselt.csv","w") |
| 112 |
t=0 |
| 113 |
t_out=0 |
| 114 |
n=0 |
| 115 |
n_out=0 |
| 116 |
dt=None |
| 117 |
a=None |
| 118 |
dt_a=None |
| 119 |
|
| 120 |
vol=integrate(1.,Function(dom)) |
| 121 |
x=dom.getX() |
| 122 |
# |
| 123 |
# set up heat problem: |
| 124 |
# |
| 125 |
heat=TemperatureCartesian(dom,theta=THETA,useSUPG=SUPG) |
| 126 |
heat.setTolerance(TOL) |
| 127 |
|
| 128 |
fixed_T_at=whereZero(x[DIM-1])+whereZero(H-x[DIM-1]) |
| 129 |
heat.setInitialTemperature(T) |
| 130 |
print "initial Temperature range ",inf(T),sup(T) |
| 131 |
heat.setValue(rhocp=Scalar(1.,Function(dom)),k=Scalar(1.,Function(dom)),given_T_mask=fixed_T_at) |
| 132 |
# |
| 133 |
# set up velovity problem |
| 134 |
# |
| 135 |
sp=StokesProblemCartesian(dom) |
| 136 |
sp.setTolerance(TOL*extratol) |
| 137 |
sp.setToleranceReductionFactor(TOL) |
| 138 |
x2=ReducedSolution(dom).getX() |
| 139 |
p=-RA*(x2[DIM-1]-0.5*x2[DIM-1]**2) |
| 140 |
|
| 141 |
fixed_v_mask=Vector(0,Solution(dom)) |
| 142 |
for d in range(DIM): |
| 143 |
if d == DIM-1: |
| 144 |
ll = H |
| 145 |
else: |
| 146 |
ll = L |
| 147 |
fixed_v_mask+=(whereZero(x[d])+whereZero(x[d]-ll))*unitVector(d,DIM) |
| 148 |
|
| 149 |
# |
| 150 |
# let the show begin: |
| 151 |
# |
| 152 |
while t<T_END: |
| 153 |
v_last=v*1 |
| 154 |
print "============== solve for v ========================" |
| 155 |
viscosity=exp(A*(1./(1+T.interpolate(Function(dom)))-1./2.)) |
| 156 |
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) |
| 158 |
#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') |
| 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): |
| 164 |
print "range %d-velocity"%d,inf(v[d]),sup(v[d]) |
| 165 |
|
| 166 |
if t>=t_out: |
| 167 |
saveVTK("state.%d.vtu"%n_out,T=T,v=v) |
| 168 |
print "visualization file %d for time step %e generated."%(n_out,t) |
| 169 |
n_out+=1 |
| 170 |
t_out+=DT_OUT |
| 171 |
Nu=1.+integrate(viscosity*length(grad(v))**2)/(RA*vol) |
| 172 |
if dom.getMPIRank() ==0: nusselt_file.write("%e %e\n"%(t,Nu)) |
| 173 |
heat.setValue(v=v,Q=DI/RA*viscosity*length(symmetric(grad(v)))**2) |
| 174 |
print "nusselt number = ",Nu,n |
| 175 |
if n>0: |
| 176 |
a,a_alt = (v_last-v)/dt, a |
| 177 |
dt_a,dt_a_alt = dt, dt_a |
| 178 |
if n>1: |
| 179 |
z=(a-a_alt)/((dt_a+dt_a_alt)/2) |
| 180 |
f=Lsup(z)/Lsup(v) |
| 181 |
print "estimated error ",f*dt**2 |
| 182 |
dt_new=min(2*dt,max(dt/2,sqrt(0.05/f))) |
| 183 |
dt=min(dt_new,heat.getSafeTimeStepSize()) |
| 184 |
else: |
| 185 |
dt=heat.getSafeTimeStepSize() |
| 186 |
print n,". time step t=",t," step size ",dt |
| 187 |
print "============== solve for T ========================" |
| 188 |
T=heat.solve(dt, verbose=VERBOSE) |
| 189 |
print "Temperature range ",inf(T),sup(T) |
| 190 |
n+=1 |
| 191 |
t+=dt |
| 192 |
# ========================= |
| 193 |
# |
| 194 |
# create restart files: |
| 195 |
# |
| 196 |
if create_restartfiles_every_step>0: |
| 197 |
if (n-1)%create_restartfiles_every_step == 0: |
| 198 |
c=(n-1)/create_restartfiles_every_step |
| 199 |
old_restart_dir="restart_%s_"%(c-1) |
| 200 |
new_restart_dir="restart_%s_"%c |
| 201 |
|
| 202 |
print "Write restart files to ",new_restart_dir |
| 203 |
if not os.path.isdir(new_restart_dir): os.mkdir(new_restart_dir) |
| 204 |
v.dump(os.path.join(new_restart_dir,"v.nc")) |
| 205 |
p.dump(os.path.join(new_restart_dir,"p.nc")) |
| 206 |
T.dump(os.path.join(new_restart_dir,"T.nc")) |
| 207 |
if n>1: |
| 208 |
file(os.path.join(new_restart_dir,"stamp.%d"%dom.getMPIRank()),"w").write("%e; %e; %s; %s; %e; %e;\n"%(t, t_out, n, n_out, dt, dt_a)) |
| 209 |
a.dump(os.path.join(new_restart_dir,"a.nc")) |
| 210 |
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)) |
| 212 |
removeRestartDirectory(old_restart_dir) |
| 213 |
elapsed = time.time() - t1 |
| 214 |
print "plot","\t",NE,"\t",elapsed |
| 215 |
|
Parent Directory
| 
Revision Log