| 1 |
from esys.escript import * |
| 2 |
from esys.escript.models import TemperatureCartesian, StokesProblemCartesian |
| 3 |
from esys.finley import Rectangle, Brick |
| 4 |
from math import pi, ceil |
| 5 |
NE=40 |
| 6 |
DIM=3 |
| 7 |
H=1. |
| 8 |
L=4*H |
| 9 |
THETA=0.5 |
| 10 |
TOL=1.e-3 |
| 11 |
PERTURBATION=0.1 |
| 12 |
T_END=0.3 |
| 13 |
DT_OUT=T_END/500 |
| 14 |
VERBOSE=False |
| 15 |
RA=1.e5 # Rayleigh number |
| 16 |
A=22. # Arenious number |
| 17 |
DI = 0. # dissipation number |
| 18 |
SUPG=False |
| 19 |
print "total number of elements = ",NE**DIM*int(L/H)**(DIM-1) |
| 20 |
# |
| 21 |
# set up domain: |
| 22 |
# |
| 23 |
if DIM==2: |
| 24 |
dom=Rectangle(int(ceil(L*NE/H)),NE,l0=L,l1=H,order=2, useFullElementOrder=True) |
| 25 |
else: |
| 26 |
dom=Brick(int(ceil(L*NE/H)),int(ceil(L*NE/H)),NE,l0=L,l1=L,l2=H,order=2, useFullElementOrder=True) |
| 27 |
|
| 28 |
vol=integrate(1.,Function(dom)) |
| 29 |
x=dom.getX() |
| 30 |
# |
| 31 |
# set up heat problem: |
| 32 |
# |
| 33 |
heat=TemperatureCartesian(dom,theta=THETA,useSUPG=SUPG) |
| 34 |
heat.setTolerance(TOL) |
| 35 |
|
| 36 |
fixed_T_at=whereZero(x[DIM-1])+whereZero(H-x[DIM-1]) |
| 37 |
T=Scalar(1,ReducedSolution(dom)) |
| 38 |
for d in range(DIM): |
| 39 |
if d == DIM-1: |
| 40 |
T*=sin(x[d]/H*pi) |
| 41 |
else: |
| 42 |
T*=cos(x[d]/H*pi) |
| 43 |
T=1.-x[DIM-1]+PERTURBATION*T |
| 44 |
heat.setInitialTemperature(T) |
| 45 |
print "initial Temperature range ",inf(T),sup(T) |
| 46 |
heat.setValue(rhocp=Scalar(1.,Function(dom)),k=Scalar(1.,Function(dom)),given_T_mask=fixed_T_at) |
| 47 |
# |
| 48 |
# set up velovity problem |
| 49 |
# |
| 50 |
sp=StokesProblemCartesian(dom) |
| 51 |
sp.setTolerance(TOL) |
| 52 |
sp.setToleranceReductionFactor(TOL) |
| 53 |
x2=ReducedSolution(dom).getX() |
| 54 |
p=-RA*(x2[DIM-1]-0.5*x2[DIM-1]**2) |
| 55 |
v=Vector(0,Solution(dom)) |
| 56 |
|
| 57 |
fixed_v_mask=Vector(0,Solution(dom)) |
| 58 |
for d in range(DIM): |
| 59 |
if d == DIM-1: |
| 60 |
ll = H |
| 61 |
else: |
| 62 |
ll = L |
| 63 |
fixed_v_mask+=(whereZero(x[d])+whereZero(x[d]-ll))*unitVector(d,DIM) |
| 64 |
|
| 65 |
# |
| 66 |
# let the show begin: |
| 67 |
# |
| 68 |
t=0 |
| 69 |
t_out=0 |
| 70 |
n=0 |
| 71 |
n_out=0 |
| 72 |
dt=None |
| 73 |
dt_new=1. |
| 74 |
a=0 |
| 75 |
if dom.getMPIRank() ==0: |
| 76 |
nusselt_file=open("nusselt.csv","w") |
| 77 |
while t<T_END: |
| 78 |
v_last=v*1. |
| 79 |
print "============== solve for v ========================" |
| 80 |
viscosity=exp(A*(1./(1+T.interpolate(Function(dom)))-0.5)) |
| 81 |
print "viscosity range :", inf(viscosity), sup(viscosity) |
| 82 |
sp.initialize(f=T*(RA*unitVector(DIM-1,DIM)),eta=viscosity,fixed_u_mask=fixed_v_mask) |
| 83 |
v,p=sp.solve(v,p,show_details=VERBOSE, verbose=True,max_iter=500) |
| 84 |
|
| 85 |
for d in range(DIM): |
| 86 |
print "range %d-velocity"%d,inf(v[d]),sup(v[d]) |
| 87 |
|
| 88 |
if t>=t_out: |
| 89 |
saveVTK("state.%d.vtu"%n_out,T=T,v=v) |
| 90 |
print "visualization file %d for time step %e generated."%(n_out,t) |
| 91 |
n_out+=1 |
| 92 |
t_out+=DT_OUT |
| 93 |
Nu=1.+integrate(viscosity*length(grad(v))**2)/(RA*vol) |
| 94 |
if dom.getMPIRank() ==0: nusselt_file.write("%e %e\n"%(t,Nu)) |
| 95 |
print "nusselt number = ",Nu |
| 96 |
if n>0: |
| 97 |
a,a_alt = (v_last-v)/dt, a |
| 98 |
dt_alt=dt |
| 99 |
if n>1: |
| 100 |
z=(a-a_alt)/((dt+dt_alt)/2) |
| 101 |
f=Lsup(z)/Lsup(v) |
| 102 |
print "estimated error ",f*dt**2 |
| 103 |
dt_new, dt_alt =min(2*dt,max(dt/2,sqrt(0.05/f))), dt |
| 104 |
# dt_new, dt_alt =sqrt(0.05/f), dt |
| 105 |
heat.setValue(v=v,Q=DI/RA*viscosity*length(symmetric(grad(v)))**2) |
| 106 |
dt=min(dt_new,heat.getSafeTimeStepSize()) |
| 107 |
print n,". time step t=",t," step size ",dt |
| 108 |
print "============== solve for T ========================" |
| 109 |
T=heat.solve(dt, verbose=VERBOSE) |
| 110 |
print "Temperature range ",inf(T),sup(T) |
| 111 |
n+=1 |
| 112 |
t+=dt |
Parent Directory
| 
Revision Log