/[ff3d]/ff3d/testsuite/natural-convection.pl
ViewVC logotype

Diff of /ff3d/testsuite/natural-convection.pl

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 1.1 by ArceneL, Fri Sep 26 07:55:59 2003 UTC revision 1.2 by ArceneL, Tue Oct 14 15:29:59 2003 UTC
# Line 24  Line 24 
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();

Legend:
Removed from v.1.1  
changed lines
  Added in v.1.2

savannah-hackers-public@gnu.org
ViewVC Help
Powered by ViewVC 1.1.26