24 |
|
|
25 |
require 'check.ph'; |
require 'check.ph'; |
26 |
|
|
27 |
# Povray scene: two cylinders assembled in a T shape. |
Run('','natural-convection.ff'); |
|
WriteFile('test.pov',<< 'EOF'); |
|
|
union { |
|
|
cylinder {<0.3,0.5,0.5>, <1,0.5,0.5>, 0.1 } |
|
|
cylinder {<0.3,0.5,0>, <0.3,0.5,1>, 0.1 } |
|
|
|
|
|
pigment { color rgb <1,0,0> } |
|
|
} |
|
|
|
|
|
EOF |
|
|
|
|
|
Run(<< 'EOF'); |
|
|
vector n = (20,20,20); |
|
|
vertex A = (0,0,0); |
|
|
vertex B = (1,1,1); |
|
|
|
|
|
scene S = pov("test.pov"); |
|
|
mesh M = structured(n,A,B); |
|
|
domain Omega = domain(S,outside(<1,0,0>)); |
|
|
mesh T=surface(Omega,M); save(opendx,"T.dat",T); |
|
|
|
|
|
function u1 = 0; |
|
|
function u2 = 0; |
|
|
function u3 = 0; |
|
|
function p = 0; |
|
|
|
|
|
function teta0 = 0; |
|
|
|
|
|
double dt = 1; |
|
|
double K = 0.1; |
|
|
double eps = 1E-3; |
|
|
|
|
|
function Ki = 1-one(<1,0,0>); |
|
|
|
|
|
double i=0; |
|
|
for (i=0; i<10; i=i+1) { |
|
|
solve(teta) in Omega by M method(type=penalty), |
|
|
krylov(precond=diagonal), memory(matrix=sparse), cg(epsilon=1E-10) |
|
|
{ |
|
|
pde(teta) |
|
|
teta-div(dt*grad(teta))=Ki*convect(Ki*u1,Ki*u2,Ki*u3,dt,teta0); |
|
|
teta = 100 on T; |
|
|
teta = 0 on M; |
|
|
}; |
|
|
|
|
|
teta0 = Ki*teta; |
|
|
u1 = 0; |
|
|
u2 = 0; |
|
|
u3 = 0; |
|
|
p = 0; |
|
|
|
|
|
solve(u1,u2,u3,p) in Omega by M |
|
|
method(type=eliminate), krylov(precond=none,type=cg), |
|
|
cg(epsilon=1E-5), memory(matrix=none) |
|
|
{ |
|
|
pde(u1) |
|
|
(1-Ki)*u1-div(Ki*grad(u1))+Ki*dx(p)=0; |
|
|
|
|
|
u1 = 0 on M; |
|
|
u1 = 0 on T; |
|
|
|
|
|
pde(u2) |
|
|
(1-Ki)*u2-div(Ki*grad(u2))+Ki*dy(p)=K*teta0; |
|
|
|
|
|
u2 = 0 on M; |
|
|
u2 = 0 on T; |
|
|
|
|
|
pde(u3) |
|
|
(1-Ki)*u3-div(Ki*grad(u3))+Ki*dz(p)=0; |
|
|
|
|
|
u3 = 0 on M; |
|
|
u3 = 0 on T; |
|
|
|
|
|
pde(p) |
|
|
-div(eps*grad(p))+eps*p+Ki*dx(u1)+Ki*dy(u2)+Ki*dz(u3)=0; |
|
|
}; |
|
|
save(opendx,"u".i.".dat", u1, M); save(opendx,"v".i.".dat", u2, M); |
|
|
save(opendx,"w".i.".dat", u3, M); save(opendx,"p".i.".dat", p, M); |
|
|
|
|
|
save(opendx,"teta".i.".dat", teta0, M); |
|
|
} |
|
|
EOF |
|
28 |
|
|
29 |
# Compares the resulting output with a reference |
# Compares the resulting output with a reference |
30 |
|
|
31 |
VectorCheck('test.bb','natural-convection-ref.bb',1e-6); |
VectorCheck('natural-convection-u-10.bb','natural-convection-u-ref.bb',1e-6); |
32 |
|
VectorCheck('natural-convection-v-10.bb','natural-convection-v-ref.bb',1e-6); |
33 |
|
VectorCheck('natural-convection-w-10.bb','natural-convection-w-ref.bb',1e-6); |
34 |
|
VectorCheck('natural-convection-p-10.bb','natural-convection-p-ref.bb',1e-6); |
35 |
|
VectorCheck('natural-convection-theta-10.bb', |
36 |
|
'natural-convection-theta-ref.bb',1e-6); |
37 |
|
|
38 |
Report(); |
Report(); |