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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 944 - (show annotations)
Tue Jan 30 08:57:37 2007 UTC (12 years, 4 months ago) by gross
File MIME type: text/x-python
File size: 4631 byte(s)
PropertySet added
1 # $Id$
2
3 __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 import sys
9 import os
10 import unittest
11 import time
12
13 from esys.escript import *
14 from esys.escript.linearPDEs import *
15 from esys import finley
16
17 starttime = time.time()
18
19 print "\nSimpleSolve.py"
20 print "--------------"
21
22 alpha=0.7
23 error_tol=1.e-5
24
25 # generate mesh
26
27 # print "\nGenerate mesh: finley.Rectangle(9,12,1)=>"
28 # mydomain=finley.Rectangle(140,140)
29
30 # print "\nGenerate mesh: finley.Rectangle(4,4,1)=>"
31 # mydomain=finley.Rectangle(5,5,1)
32 mydomain=finley.Rectangle(500,500,1)
33 mydomain=finley.Rectangle(150,150,1)
34
35 print "\nGenerate mesh: finley.Rectangle(151,151,1)=>"
36 # mydomain=finley.Rectangle(151,151,1)
37 # mydomain=finley.Rectangle(128,128,1)
38
39 # set the direct solver switch
40 DIRECT=LinearPDE.DIRECT
41 # DIRECT=LinearPDE.ITERATIVE # this will switch of the DIRECT solver call to avoid external library calls which may not be available everywhere
42
43 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 print "is x protected: ",x.isProtected()
63
64
65 print "norm_u_ex=Lsup(u_ex):"
66 norm_u_ex=Lsup(u_ex)
67 print "is u_ex protected: ",u_ex.isProtected()
68
69 print "\nGenerate a test solution (1)"
70 print "----------------------------"
71
72 print "mypde=LinearPDE( A=[[1.,0.8],[0.4,1.]], D=alpha, Y=alpha, domain=mydomain)"
73 mypde=LinearPDE(mydomain)
74 mypde.setDebugOn()
75 mypde.setValue(A=[[1.,-0.001],[-0.001,1.]],D=alpha,Y=alpha)
76
77 print "mypde.checkSymmetry()"
78 print mypde.checkSymmetry()
79
80 print "\nIterative Solver (1)=>"
81 # mypde.setSolverMethod(mypde.PRES20,preconditioner=mypde.ILU0)
82 mypde.setSolverMethod(mypde.BICGSTAB,preconditioner=mypde.JACOBI)
83 u_i=mypde.getSolution(verbose=True,iter_max=3000)
84
85 print "\nDirect Solver (1)=>"
86 mypde.setSolverMethod(DIRECT)
87 u_d=mypde.getSolution(verbose=True)
88
89 print "\n***************************************************************"
90 error=u_ex-u_d
91 error_norm=Lsup(error)/norm_u_ex
92 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 error=u_ex-u_i
97 error_norm=Lsup(error)/norm_u_ex
98 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 print "***************************************************************"
103 # del mypde
104 del u_i
105 # releaseUnusedMemory()
106 print "***************************************************************"
107
108
109
110 # 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 print " msk=whereZero(x[0])+whereZero(x[0]-1.)"
119 msk=whereZero(x[0])+whereZero(x[0]-1.)
120
121 print "mypde=LinearPDE(A=[[1.,0.],[0.,1.]],q=msk,r=u_ex)"
122 mypde=LinearPDE(mydomain)
123 mypde.setDebugOn()
124 mypde.setValue(A=[[1.,0.0],[0.0,1.]],q=msk,r=u_ex)
125
126 print "mypde.checkSymmetry()"
127 print mypde.checkSymmetry()
128
129 # generate a test solution 2
130
131 print "\nGenerate a test solution (2)"
132 print "----------------------------"
133
134 print "\nDirect Solver (2)=>"
135
136 mypde.setSymmetryOn()
137 mypde.setTolerance(1.e-13)
138
139 # mypde.setSymmetryOn() : is not woking yet!
140 mypde.setSolverMethod(DIRECT)
141 u_d=mypde.getSolution(verbose=True)
142
143 print "\nIterative Solver (2)=>"
144
145 mypde.setSolverMethod(mypde.PCG)
146 u_i=mypde.getSolution(verbose=True,iter_max=3000)
147
148 print "\n******************************************************************"
149 error=u_ex-u_d
150 error_norm=Lsup(error)/norm_u_ex
151 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 error=u_ex-u_i
156 error_norm=Lsup(error)/norm_u_ex
157 print "norm of the error for iterative solver is: ",error_norm
158 if error_norm > error_tol:
159 print "### error norm exceeded maximum tolerance ###"
160 sys.exit(1)
161 print "******************************************************************"
162
163 print "\n-----"
164 print "Done."
165 print "-----"
166
167 stoptime = time.time()
168 elapsed = stoptime - starttime
169 print "\nElapsed time: ", elapsed, "\n"
170
171 sys.exit(0)

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26