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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 944 - (hide annotations)
Tue Jan 30 08:57:37 2007 UTC (13 years, 2 months ago) by gross
File MIME type: text/x-python
File size: 4631 byte(s)
PropertySet added
1 jgs 82 # $Id$
2    
3 elspeth 617 __copyright__=""" Copyright (c) 2006 by ACcESS MNRF
4     http://www.access.edu.au
5     Primary Business: Queensland, Australia"""
6     __license__="""Licensed under the Open Software License version 3.0
7     http://www.opensource.org/licenses/osl-3.0.php"""
8 jgs 82 import sys
9     import os
10     import unittest
11 jgs 153 import time
12 jgs 82
13 jgs 149 from esys.escript import *
14     from esys.escript.linearPDEs import *
15     from esys import finley
16 jgs 82
17 gross 853 starttime = time.time()
18 jgs 153
19 jgs 82 print "\nSimpleSolve.py"
20     print "--------------"
21    
22 jgs 150 alpha=0.7
23 jgs 148 error_tol=1.e-5
24 jgs 82
25     # generate mesh
26    
27 jgs 102 # print "\nGenerate mesh: finley.Rectangle(9,12,1)=>"
28 jgs 115 # mydomain=finley.Rectangle(140,140)
29 jgs 82
30 jgs 102 # print "\nGenerate mesh: finley.Rectangle(4,4,1)=>"
31 gross 806 # mydomain=finley.Rectangle(5,5,1)
32 gross 777 mydomain=finley.Rectangle(500,500,1)
33 gross 944 mydomain=finley.Rectangle(150,150,1)
34 jgs 82
35     print "\nGenerate mesh: finley.Rectangle(151,151,1)=>"
36 jgs 148 # mydomain=finley.Rectangle(151,151,1)
37 jgs 115 # mydomain=finley.Rectangle(128,128,1)
38 jgs 82
39 gross 778 # set the direct solver switch
40     DIRECT=LinearPDE.DIRECT
41 gross 806 # DIRECT=LinearPDE.ITERATIVE # this will switch of the DIRECT solver call to avoid external library calls which may not be available everywhere
42 gross 778
43 jgs 82 print "\nSetup domain and functions"
44     print "--------------------------"
45    
46     print "e=Function(mydomain):"
47     e=Function(mydomain)
48    
49     print "n=ContinuousFunction(mydomain):"
50     n=ContinuousFunction(mydomain)
51    
52     # get handles to nodes and elements 1
53    
54     print "\nGet handles to nodes and elements(1)=>"
55     print "--------------------------------------"
56    
57     print "u_ex=Scalar(1,n,True):"
58     u_ex=Scalar(1,n,True)
59    
60     print "x=e.getX():"
61     x=e.getX()
62 gross 783 print "is x protected: ",x.isProtected()
63 jgs 82
64 gross 783
65 gross 298 print "norm_u_ex=Lsup(u_ex):"
66     norm_u_ex=Lsup(u_ex)
67 gross 783 print "is u_ex protected: ",u_ex.isProtected()
68 jgs 82
69 jgs 108 print "\nGenerate a test solution (1)"
70     print "----------------------------"
71    
72 jgs 102 print "mypde=LinearPDE( A=[[1.,0.8],[0.4,1.]], D=alpha, Y=alpha, domain=mydomain)"
73 jgs 104 mypde=LinearPDE(mydomain)
74 jgs 102 mypde.setDebugOn()
75 gross 432 mypde.setValue(A=[[1.,-0.001],[-0.001,1.]],D=alpha,Y=alpha)
76 jgs 82
77 jgs 108 print "mypde.checkSymmetry()"
78     print mypde.checkSymmetry()
79 jgs 82
80 jgs 102 print "\nIterative Solver (1)=>"
81 gross 453 # mypde.setSolverMethod(mypde.PRES20,preconditioner=mypde.ILU0)
82     mypde.setSolverMethod(mypde.BICGSTAB,preconditioner=mypde.JACOBI)
83 gross 431 u_i=mypde.getSolution(verbose=True,iter_max=3000)
84 jgs 102
85 jgs 100 print "\nDirect Solver (1)=>"
86 gross 778 mypde.setSolverMethod(DIRECT)
87 jgs 149 u_d=mypde.getSolution(verbose=True)
88 jgs 97
89 jgs 82 print "\n***************************************************************"
90     error=u_ex-u_d
91 gross 298 error_norm=Lsup(error)/norm_u_ex
92 jgs 147 print "norm of the error for direct solver is : ",error_norm
93     if error_norm > error_tol:
94     print "### error norm exceeded maximum tolerance ###"
95     sys.exit(1)
96 jgs 82 error=u_ex-u_i
97 gross 298 error_norm=Lsup(error)/norm_u_ex
98 jgs 147 print "norm of the error for iterative solver is: ",error_norm
99     if error_norm > error_tol:
100     print "### error norm exceeded maximum tolerance ###"
101     sys.exit(1)
102 jgs 82 print "***************************************************************"
103 gross 820 # del mypde
104     del u_i
105     # releaseUnusedMemory()
106 gross 425 print "***************************************************************"
107 jgs 82
108 gross 425
109 gross 820
110 jgs 82 # get handles to nodes and elements 2
111    
112     print "\nGet handles to nodes and elements(2)=>"
113     print "--------------------------------------"
114    
115     print "x=n.getX():"
116     x=n.getX()
117    
118 gross 298 print " msk=whereZero(x[0])+whereZero(x[0]-1.)"
119     msk=whereZero(x[0])+whereZero(x[0]-1.)
120 jgs 82
121 jgs 102 print "mypde=LinearPDE(A=[[1.,0.],[0.,1.]],q=msk,r=u_ex)"
122 jgs 104 mypde=LinearPDE(mydomain)
123 jgs 108 mypde.setDebugOn()
124 gross 432 mypde.setValue(A=[[1.,0.0],[0.0,1.]],q=msk,r=u_ex)
125 jgs 82
126 jgs 108 print "mypde.checkSymmetry()"
127     print mypde.checkSymmetry()
128    
129 jgs 82 # generate a test solution 2
130    
131     print "\nGenerate a test solution (2)"
132     print "----------------------------"
133    
134     print "\nDirect Solver (2)=>"
135    
136 gross 404 mypde.setSymmetryOn()
137 jgs 123 mypde.setTolerance(1.e-13)
138 jgs 115
139 jgs 102 # mypde.setSymmetryOn() : is not woking yet!
140 gross 778 mypde.setSolverMethod(DIRECT)
141 jgs 149 u_d=mypde.getSolution(verbose=True)
142 jgs 82
143     print "\nIterative Solver (2)=>"
144    
145 gross 778 mypde.setSolverMethod(mypde.PCG)
146 jgs 149 u_i=mypde.getSolution(verbose=True,iter_max=3000)
147 jgs 82
148     print "\n******************************************************************"
149     error=u_ex-u_d
150 gross 298 error_norm=Lsup(error)/norm_u_ex
151 jgs 147 print "norm of the error for direct solver is : ",error_norm
152     if error_norm > error_tol:
153     print "### error norm exceeded maximum tolerance ###"
154     sys.exit(1)
155 jgs 82 error=u_ex-u_i
156 gross 298 error_norm=Lsup(error)/norm_u_ex
157 jgs 147 print "norm of the error for iterative solver is: ",error_norm
158 jgs 148 if error_norm > error_tol:
159 jgs 147 print "### error norm exceeded maximum tolerance ###"
160     sys.exit(1)
161 jgs 82 print "******************************************************************"
162    
163     print "\n-----"
164     print "Done."
165     print "-----"
166 jgs 147
167 gross 853 stoptime = time.time()
168 jgs 153 elapsed = stoptime - starttime
169     print "\nElapsed time: ", elapsed, "\n"
170    
171 jgs 147 sys.exit(0)

Properties

Name Value
svn:eol-style native
svn:keywords Author Date Id Revision

  ViewVC Help
Powered by ViewVC 1.1.26