/[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 778 - (show annotations)
Thu Jul 13 02:27:45 2006 UTC (13 years ago) by gross
File MIME type: text/x-python
File size: 4513 byte(s)
add a switch to stop calls of external solver libraries
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.clock()
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(50,50,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
63 print "norm_u_ex=Lsup(u_ex):"
64 norm_u_ex=Lsup(u_ex)
65
66 print "\nGenerate a test solution (1)"
67 print "----------------------------"
68
69 print "mypde=LinearPDE( A=[[1.,0.8],[0.4,1.]], D=alpha, Y=alpha, domain=mydomain)"
70 mypde=LinearPDE(mydomain)
71 mypde.setDebugOn()
72 mypde.setValue(A=[[1.,-0.001],[-0.001,1.]],D=alpha,Y=alpha)
73
74 print "mypde.checkSymmetry()"
75 print mypde.checkSymmetry()
76
77 print "\nIterative Solver (1)=>"
78 # mypde.setSolverMethod(mypde.PRES20,preconditioner=mypde.ILU0)
79 mypde.setSolverMethod(mypde.BICGSTAB,preconditioner=mypde.JACOBI)
80 u_i=mypde.getSolution(verbose=True,iter_max=3000)
81
82 print "\nDirect Solver (1)=>"
83 mypde.setSolverMethod(DIRECT)
84 u_d=mypde.getSolution(verbose=True)
85
86 print "\n***************************************************************"
87 error=u_ex-u_d
88 error_norm=Lsup(error)/norm_u_ex
89 print "norm of the error for direct solver is : ",error_norm
90 if error_norm > error_tol:
91 print "### error norm exceeded maximum tolerance ###"
92 sys.exit(1)
93 error=u_ex-u_i
94 error_norm=Lsup(error)/norm_u_ex
95 print "norm of the error for iterative solver is: ",error_norm
96 if error_norm > error_tol:
97 print "### error norm exceeded maximum tolerance ###"
98 sys.exit(1)
99 print "***************************************************************"
100 del mypde
101 print "***************************************************************"
102
103
104 # get handles to nodes and elements 2
105
106 print "\nGet handles to nodes and elements(2)=>"
107 print "--------------------------------------"
108
109 print "x=n.getX():"
110 x=n.getX()
111
112 print " msk=whereZero(x[0])+whereZero(x[0]-1.)"
113 msk=whereZero(x[0])+whereZero(x[0]-1.)
114
115 print "mypde=LinearPDE(A=[[1.,0.],[0.,1.]],q=msk,r=u_ex)"
116 mypde=LinearPDE(mydomain)
117 mypde.setDebugOn()
118 mypde.setValue(A=[[1.,0.0],[0.0,1.]],q=msk,r=u_ex)
119
120 print "mypde.checkSymmetry()"
121 print mypde.checkSymmetry()
122
123 # generate a test solution 2
124
125 print "\nGenerate a test solution (2)"
126 print "----------------------------"
127
128 print "\nDirect Solver (2)=>"
129
130 mypde.setSymmetryOn()
131 mypde.setTolerance(1.e-13)
132
133 # mypde.setSymmetryOn() : is not woking yet!
134 mypde.setSolverMethod(DIRECT)
135 u_d=mypde.getSolution(verbose=True)
136
137 print "\nIterative Solver (2)=>"
138
139 mypde.setSolverMethod(mypde.PCG)
140 u_i=mypde.getSolution(verbose=True,iter_max=3000)
141
142 print "\n******************************************************************"
143 error=u_ex-u_d
144 error_norm=Lsup(error)/norm_u_ex
145 print "norm of the error for direct solver is : ",error_norm
146 if error_norm > error_tol:
147 print "### error norm exceeded maximum tolerance ###"
148 sys.exit(1)
149 error=u_ex-u_i
150 error_norm=Lsup(error)/norm_u_ex
151 print "norm of the error for iterative solver is: ",error_norm
152 if error_norm > error_tol:
153 print "### error norm exceeded maximum tolerance ###"
154 sys.exit(1)
155 print "******************************************************************"
156
157 print "\n-----"
158 print "Done."
159 print "-----"
160
161 stoptime = time.clock()
162 elapsed = stoptime - starttime
163 print "\nElapsed time: ", elapsed, "\n"
164
165 sys.exit(0)

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26