/[escript]/trunk/finley/test/python/SimpleSolve.py
ViewVC logotype

Diff of /trunk/finley/test/python/SimpleSolve.py

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

revision 82 by jgs, Tue Oct 26 06:53:54 2004 UTC revision 148 by jgs, Tue Aug 23 01:24:31 2005 UTC
# Line 4  import sys Line 4  import sys
4  import os  import os
5  import unittest  import unittest
6    
7  esys_root=os.getenv('ESYS_ROOT')  from escript.escript import *
8  sys.path.append(esys_root+'/finley/lib')  from escript.linearPDEs import *
9  sys.path.append(esys_root+'/escript/lib')  from finley import finley
 sys.path.append(esys_root+'/escript/py_src')  
   
 from escript import *  
 from util import *  
 from linearPDE import *  
 import finley  
10    
11  print "\nSimpleSolve.py"  print "\nSimpleSolve.py"
12  print "--------------"  print "--------------"
13    
14  my_options = {  alpha=0.07
15            "verbose" : True,  error_tol=1.e-5
           "reordering" : NO_REORDERING,  
           "tolerance" : 1.E-8,  
           "final_residual" : 0.,  
           "iterative_method" : PCG,  
           "preconditioner" : JACOBI,  
           "iter_max" :  5000,  
           "iter" : 0,  
           "drop_tolerance" : 1.10,  
           "drop_storage" : 2.  
 }  
   
 alpha=0.01  
   
 print "\nOptions: ", my_options  
16    
17  # generate mesh  # generate mesh
18    
19  print "\nGenerate mesh: finley.Rectangle(9,12,1)=>"  # print "\nGenerate mesh: finley.Rectangle(9,12,1)=>"
20  mydomain=finley.Rectangle(9,12,1)  # mydomain=finley.Rectangle(140,140)
21    
22  print "\nGenerate mesh: finley.Rectangle(4,4,1)=>"  # print "\nGenerate mesh: finley.Rectangle(4,4,1)=>"
23  mydomain=finley.Rectangle(4,4,1)  mydomain=finley.Rectangle(190,190,1)
24    
25  print "\nGenerate mesh: finley.Rectangle(151,151,1)=>"  print "\nGenerate mesh: finley.Rectangle(151,151,1)=>"
26  mydomain=finley.Rectangle(151,151,1)  # mydomain=finley.Rectangle(151,151,1)
27    # mydomain=finley.Rectangle(128,128,1)
28    
29  print "\nSetup domain and functions"  print "\nSetup domain and functions"
30  print "--------------------------"  print "--------------------------"
# Line 68  x=e.getX() Line 49  x=e.getX()
49  print "norm_u_ex=u_ex.Lsup():"  print "norm_u_ex=u_ex.Lsup():"
50  norm_u_ex=u_ex.Lsup()  norm_u_ex=u_ex.Lsup()
51    
 print "mypde=linearPDE( A=[[1.,0.7],[0.7,1.]], D=alpha, Y=alpha, domain=mydomain)"  
 mypde=linearPDE(A=[[1.,0.7],[0.7,1.]],D=alpha,Y=alpha,domain=mydomain)  
 #mypde=linearPDE(A=[[1.,0.],[0.,1.]],D=alpha,Y=alpha,domain=mydomain)  
   
 # generate a test solution 1  
   
52  print "\nGenerate a test solution (1)"  print "\nGenerate a test solution (1)"
53  print "----------------------------"  print "----------------------------"
54    
55  print "\nDirect Solver (1)=>"  print "mypde=LinearPDE( A=[[1.,0.8],[0.4,1.]], D=alpha, Y=alpha, domain=mydomain)"
56    mypde=LinearPDE(mydomain)
57    mypde.setDebugOn()
58    mypde.setValue(A=[[1.,0.8],[0.4,1.]],D=alpha,Y=alpha)
59    
60  u_d=mypde.getSolution(iterative=False,**my_options)  print "mypde.checkSymmetry()"
61    print mypde.checkSymmetry()
62    
63  print "\nIterative Solver (1)=>"  print "\nIterative Solver (1)=>"
64    # u_i=mypde.getSolution(preconditioner=ILU0,iter_max=3000)
65    u_i=mypde.getSolution(iter_max=3000)
66    
67  u_i=mypde.getSolution(iterative=True,**my_options)  print "\nDirect Solver (1)=>"
68    mypde.setSolverMethod(mypde.DIRECT)
69    u_d=mypde.getSolution()
70    
71  print "\n***************************************************************"  print "\n***************************************************************"
72  error=u_ex-u_d  error=u_ex-u_d
73  print "norm of the error for direct solver is   : ",error.Lsup()/norm_u_ex  error_norm=error.Lsup()/norm_u_ex
74    print "norm of the error for direct solver is   : ",error_norm
75    if error_norm > error_tol:
76      print "### error norm exceeded maximum tolerance ###"
77      sys.exit(1)
78  error=u_ex-u_i  error=u_ex-u_i
79  print "norm of the error for iterative solver is: ",error.Lsup()/norm_u_ex  error_norm=error.Lsup()/norm_u_ex
80    print "norm of the error for iterative solver is: ",error_norm
81    if error_norm > error_tol:
82      print "### error norm exceeded maximum tolerance ###"
83      sys.exit(1)
84  print "***************************************************************"  print "***************************************************************"
85    
86  # get handles to nodes and elements 2  # get handles to nodes and elements 2
# Line 103  x=n.getX() Line 94  x=n.getX()
94  print "msk=x[0].whereZero()+(x[0]-1.).whereZero()"  print "msk=x[0].whereZero()+(x[0]-1.).whereZero()"
95  msk=x[0].whereZero()+(x[0]-1.).whereZero()  msk=x[0].whereZero()+(x[0]-1.).whereZero()
96    
97  print "mypde=linearPDE(A=[[1.,0.],[0.,1.]],q=msk,r=u_ex)"  print "mypde=LinearPDE(A=[[1.,0.],[0.,1.]],q=msk,r=u_ex)"
98  mypde=linearPDE(A=[[1.,0.],[0.,1.]],q=msk,r=u_ex)  mypde=LinearPDE(mydomain)
99    mypde.setDebugOn()
100    mypde.setValue(A=[[1.,0.],[0.,1.]],q=msk,r=u_ex)
101    
102    print "mypde.checkSymmetry()"
103    print mypde.checkSymmetry()
104    
105  # generate a test solution 2  # generate a test solution 2
106    
# Line 113  print "----------------------------" Line 109  print "----------------------------"
109    
110  print "\nDirect Solver (2)=>"  print "\nDirect Solver (2)=>"
111    
112  u_d=mypde.getSolution(iterative=False,**my_options)  #mypde.setSymmetryOn()
113    mypde.setTolerance(1.e-13)
114    
115    # mypde.setSymmetryOn() : is not woking yet!
116    mypde.setSolverMethod(mypde.DIRECT)
117    u_d=mypde.getSolution()
118    
119  print "\nIterative Solver (2)=>"  print "\nIterative Solver (2)=>"
120    
121  u_i=mypde.getSolution(iterative=True,**my_options)  mypde.setSymmetryOn()
122    mypde.setSolverMethod(mypde.DEFAULT_METHOD)
123    u_i=mypde.getSolution(iter_max=3000)
124    
125  print "\n******************************************************************"  print "\n******************************************************************"
126  error=u_ex-u_d  error=u_ex-u_d
127  print "norm of the error for direct solver is   : ",error.Lsup()/norm_u_ex  error_norm=error.Lsup()/norm_u_ex
128    print "norm of the error for direct solver is   : ",error_norm
129    if error_norm > error_tol:
130      print "### error norm exceeded maximum tolerance ###"
131      sys.exit(1)
132  error=u_ex-u_i  error=u_ex-u_i
133  print "norm of the error for iterative solver is: ",error.Lsup()/norm_u_ex  error_norm=error.Lsup()/norm_u_ex
134    print "norm of the error for iterative solver is: ",error_norm
135    if error_norm >  error_tol:
136      print "### error norm exceeded maximum tolerance ###"
137      sys.exit(1)
138  print "******************************************************************"  print "******************************************************************"
139    
140  print "\n-----"  print "\n-----"
141  print "Done."  print "Done."
142  print "-----"  print "-----"
143    
144    sys.exit(0)

Legend:
Removed from v.82  
changed lines
  Added in v.148

  ViewVC Help
Powered by ViewVC 1.1.26