/[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 617 - (show annotations)
Wed Mar 22 02:58:17 2006 UTC (14 years, 9 months ago) by elspeth
File MIME type: text/x-python
File size: 4330 byte(s)
More copyright.

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 print "\nSetup domain and functions"
40 print "--------------------------"
41
42 print "e=Function(mydomain):"
43 e=Function(mydomain)
44
45 print "n=ContinuousFunction(mydomain):"
46 n=ContinuousFunction(mydomain)
47
48 # get handles to nodes and elements 1
49
50 print "\nGet handles to nodes and elements(1)=>"
51 print "--------------------------------------"
52
53 print "u_ex=Scalar(1,n,True):"
54 u_ex=Scalar(1,n,True)
55
56 print "x=e.getX():"
57 x=e.getX()
58
59 print "norm_u_ex=Lsup(u_ex):"
60 norm_u_ex=Lsup(u_ex)
61
62 print "\nGenerate a test solution (1)"
63 print "----------------------------"
64
65 print "mypde=LinearPDE( A=[[1.,0.8],[0.4,1.]], D=alpha, Y=alpha, domain=mydomain)"
66 mypde=LinearPDE(mydomain)
67 mypde.setDebugOn()
68 mypde.setValue(A=[[1.,-0.001],[-0.001,1.]],D=alpha,Y=alpha)
69
70 print "mypde.checkSymmetry()"
71 print mypde.checkSymmetry()
72
73 print "\nIterative Solver (1)=>"
74 # mypde.setSolverMethod(mypde.PRES20,preconditioner=mypde.ILU0)
75 mypde.setSolverMethod(mypde.BICGSTAB,preconditioner=mypde.JACOBI)
76 u_i=mypde.getSolution(verbose=True,iter_max=3000)
77
78 print "\nDirect Solver (1)=>"
79 mypde.setSolverMethod(mypde.DIRECT)
80 u_d=mypde.getSolution(verbose=True)
81
82 print "\n***************************************************************"
83 error=u_ex-u_d
84 error_norm=Lsup(error)/norm_u_ex
85 print "norm of the error for direct solver is : ",error_norm
86 if error_norm > error_tol:
87 print "### error norm exceeded maximum tolerance ###"
88 sys.exit(1)
89 error=u_ex-u_i
90 error_norm=Lsup(error)/norm_u_ex
91 print "norm of the error for iterative solver is: ",error_norm
92 if error_norm > error_tol:
93 print "### error norm exceeded maximum tolerance ###"
94 sys.exit(1)
95 print "***************************************************************"
96 del mypde
97 print "***************************************************************"
98
99
100 # get handles to nodes and elements 2
101
102 print "\nGet handles to nodes and elements(2)=>"
103 print "--------------------------------------"
104
105 print "x=n.getX():"
106 x=n.getX()
107
108 print " msk=whereZero(x[0])+whereZero(x[0]-1.)"
109 msk=whereZero(x[0])+whereZero(x[0]-1.)
110
111 print "mypde=LinearPDE(A=[[1.,0.],[0.,1.]],q=msk,r=u_ex)"
112 mypde=LinearPDE(mydomain)
113 mypde.setDebugOn()
114 mypde.setValue(A=[[1.,0.0],[0.0,1.]],q=msk,r=u_ex)
115
116 print "mypde.checkSymmetry()"
117 print mypde.checkSymmetry()
118
119 # generate a test solution 2
120
121 print "\nGenerate a test solution (2)"
122 print "----------------------------"
123
124 print "\nDirect Solver (2)=>"
125
126 mypde.setSymmetryOn()
127 mypde.setTolerance(1.e-13)
128
129 # mypde.setSymmetryOn() : is not woking yet!
130 mypde.setSolverMethod(mypde.DIRECT)
131 u_d=mypde.getSolution(verbose=True)
132
133 print "\nIterative Solver (2)=>"
134
135 mypde.setSolverMethod(mypde.ITERATIVE)
136 u_i=mypde.getSolution(verbose=True,iter_max=3000)
137
138 print "\n******************************************************************"
139 error=u_ex-u_d
140 error_norm=Lsup(error)/norm_u_ex
141 print "norm of the error for direct solver is : ",error_norm
142 if error_norm > error_tol:
143 print "### error norm exceeded maximum tolerance ###"
144 sys.exit(1)
145 error=u_ex-u_i
146 error_norm=Lsup(error)/norm_u_ex
147 print "norm of the error for iterative solver is: ",error_norm
148 if error_norm > error_tol:
149 print "### error norm exceeded maximum tolerance ###"
150 sys.exit(1)
151 print "******************************************************************"
152
153 print "\n-----"
154 print "Done."
155 print "-----"
156
157 stoptime = time.clock()
158 elapsed = stoptime - starttime
159 print "\nElapsed time: ", elapsed, "\n"
160
161 sys.exit(0)

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26