/[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

trunk/esys2/finley/test/python/SimpleSolve.py revision 97 by jgs, Tue Dec 14 05:39:33 2004 UTC trunk/finley/test/python/SimpleSolve.py revision 431 by gross, Fri Jan 13 05:07:10 2006 UTC
# Line 3  Line 3 
3  import sys  import sys
4  import os  import os
5  import unittest  import unittest
6    import time
7    
8  from esys.escript import *  from esys.escript import *
9  from esys.linearPDEs import *  from esys.escript.linearPDEs import *
10  from esys import finley  from esys import finley
11    
12    starttime = time.clock()
13    
14  print "\nSimpleSolve.py"  print "\nSimpleSolve.py"
15  print "--------------"  print "--------------"
16    
17  alpha=0.01  alpha=0.7
18    error_tol=1.e-5
19    
20  # generate mesh  # generate mesh
21    
22  # print "\nGenerate mesh: finley.Rectangle(9,12,1)=>"  # print "\nGenerate mesh: finley.Rectangle(9,12,1)=>"
23  # mydomain=finley.Rectangle(1,1)  # mydomain=finley.Rectangle(140,140)
24    
25  # print "\nGenerate mesh: finley.Rectangle(4,4,1)=>"  # print "\nGenerate mesh: finley.Rectangle(4,4,1)=>"
26  # mydomain=finley.Rectangle(4,4,1)  mydomain=finley.Rectangle(150,10,1)
27    mydomain=finley.Rectangle(250,250,1)
28    # mydomain=finley.Rectangle(190,190,1)
29    
30  print "\nGenerate mesh: finley.Rectangle(151,151,1)=>"  print "\nGenerate mesh: finley.Rectangle(151,151,1)=>"
31  mydomain=finley.Rectangle(151,151,1)  # mydomain=finley.Rectangle(151,151,1)
32    # mydomain=finley.Rectangle(128,128,1)
33    
34  print "\nSetup domain and functions"  print "\nSetup domain and functions"
35  print "--------------------------"  print "--------------------------"
# Line 44  u_ex=Scalar(1,n,True) Line 51  u_ex=Scalar(1,n,True)
51  print "x=e.getX():"  print "x=e.getX():"
52  x=e.getX()  x=e.getX()
53    
54  print "norm_u_ex=u_ex.Lsup():"  print "norm_u_ex=Lsup(u_ex):"
55  norm_u_ex=u_ex.Lsup()  norm_u_ex=Lsup(u_ex)
56    
57    print "\nGenerate a test solution (1)"
58    print "----------------------------"
59    
60  print "mypde=LinearPDE( A=[[1.,0.8],[0.4,1.]], D=alpha, Y=alpha, domain=mydomain)"  print "mypde=LinearPDE( A=[[1.,0.8],[0.4,1.]], D=alpha, Y=alpha, domain=mydomain)"
61  mypde=LinearPDE(A=[[1.,0.8],[0.4,1.]],D=alpha,Y=alpha,domain=mydomain)  mypde=LinearPDE(mydomain)
62  mypde.setDebugOn()  mypde.setDebugOn()
63  #mypde=LinearPDE(A=[[1.,0.],[0.,1.]],D=alpha,Y=alpha,domain=mydomain)  mypde.setValue(A=[[1.,0.1],[0.04,1.]],D=alpha,Y=alpha)
 mypde.getOperator().saveMM("t.msh")  
64    
65  # generate a test solution 1  print "mypde.checkSymmetry()"
66    print mypde.checkSymmetry()
 print "\nGenerate a test solution (1)"  
 print "----------------------------"  
67    
68  print "\nIterative Solver (1)=>"  print "\nIterative Solver (1)=>"
69    mypde.setSolverMethod(mypde.BICGSTAB,preconditioner=mypde.ILU0)
70  u_i=mypde.getSolution()  u_i=mypde.getSolution(verbose=True,iter_max=3000)
71    
72  print "\nDirect Solver (1)=>"  print "\nDirect Solver (1)=>"
73    mypde.setSolverMethod(mypde.DIRECT)
74  mypde.setSolverMethod(DIRECT)  u_d=mypde.getSolution(verbose=True)
 u_d=mypde.getSolution()  
   
75    
76  print "\n***************************************************************"  print "\n***************************************************************"
77  error=u_ex-u_d  error=u_ex-u_d
78  print "norm of the error for direct solver is   : ",error.Lsup()/norm_u_ex  error_norm=Lsup(error)/norm_u_ex
79    print "norm of the error for direct solver is   : ",error_norm
80    if error_norm > error_tol:
81      print "### error norm exceeded maximum tolerance ###"
82      sys.exit(1)
83  error=u_ex-u_i  error=u_ex-u_i
84  print "norm of the error for iterative solver is: ",error.Lsup()/norm_u_ex  error_norm=Lsup(error)/norm_u_ex
85    print "norm of the error for iterative solver is: ",error_norm
86    if error_norm > error_tol:
87      print "### error norm exceeded maximum tolerance ###"
88      sys.exit(1)
89  print "***************************************************************"  print "***************************************************************"
90    del mypde
91    print "***************************************************************"
92    
93    
94  # get handles to nodes and elements 2  # get handles to nodes and elements 2
95    
# Line 83  print "--------------------------------- Line 99  print "---------------------------------
99  print "x=n.getX():"  print "x=n.getX():"
100  x=n.getX()  x=n.getX()
101    
102  print "msk=x[0].whereZero()+(x[0]-1.).whereZero()"  print " msk=whereZero(x[0])+whereZero(x[0]-1.)"
103  msk=x[0].whereZero()+(x[0]-1.).whereZero()  msk=whereZero(x[0])+whereZero(x[0]-1.)
104    
105  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)"
106  mypde=LinearPDE(A=[[1.,0.],[0.,1.]],q=msk,r=u_ex)  mypde=LinearPDE(mydomain)
107  mypde.setDebugOn()  mypde.setDebugOn()
108    mypde.setValue(A=[[1.,0.],[0.,1.]],q=msk,r=u_ex)
109    
110    print "mypde.checkSymmetry()"
111    print mypde.checkSymmetry()
112    
113  # generate a test solution 2  # generate a test solution 2
114    
# Line 97  print "----------------------------" Line 117  print "----------------------------"
117    
118  print "\nDirect Solver (2)=>"  print "\nDirect Solver (2)=>"
119    
120    mypde.setSymmetryOn()
121    mypde.setTolerance(1.e-13)
122    
123  # mypde.setSymmetryOn() : is not woking yet!  # mypde.setSymmetryOn() : is not woking yet!
124  mypde.setSolverMethod(DIRECT)  mypde.setSolverMethod(mypde.DIRECT)
125  u_d=mypde.getSolution()  u_d=mypde.getSolution(verbose=True)
126    
127  print "\nIterative Solver (2)=>"  print "\nIterative Solver (2)=>"
128    
129  #mypde.setSymmetryOn()  mypde.setSolverMethod(mypde.ITERATIVE)
130  mypde.setSolverMethod(DEFAULT_METHOD)  u_i=mypde.getSolution(verbose=True,iter_max=3000)
 u_i=mypde.getSolution()  
131    
132  print "\n******************************************************************"  print "\n******************************************************************"
133  error=u_ex-u_d  error=u_ex-u_d
134  print "norm of the error for direct solver is   : ",error.Lsup()/norm_u_ex  error_norm=Lsup(error)/norm_u_ex
135    print "norm of the error for direct solver is   : ",error_norm
136    if error_norm > error_tol:
137      print "### error norm exceeded maximum tolerance ###"
138      sys.exit(1)
139  error=u_ex-u_i  error=u_ex-u_i
140  print "norm of the error for iterative solver is: ",error.Lsup()/norm_u_ex  error_norm=Lsup(error)/norm_u_ex
141    print "norm of the error for iterative solver is: ",error_norm
142    if error_norm >  error_tol:
143      print "### error norm exceeded maximum tolerance ###"
144      sys.exit(1)
145  print "******************************************************************"  print "******************************************************************"
146    
147  print "\n-----"  print "\n-----"
148  print "Done."  print "Done."
149  print "-----"  print "-----"
150    
151    stoptime = time.clock()
152    elapsed = stoptime - starttime
153    print "\nElapsed time: ", elapsed, "\n"
154    
155    sys.exit(0)

Legend:
Removed from v.97  
changed lines
  Added in v.431

  ViewVC Help
Powered by ViewVC 1.1.26