/[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 404 - (show annotations)
Thu Dec 22 23:04:55 2005 UTC (13 years, 9 months ago) by gross
File MIME type: text/x-python
File size: 3849 byte(s)
direct solver works for symmetry flag now
1 # $Id$
2
3 import sys
4 import os
5 import unittest
6 import time
7
8 from esys.escript import *
9 from esys.escript.linearPDEs import *
10 from esys import finley
11
12 starttime = time.clock()
13
14 print "\nSimpleSolve.py"
15 print "--------------"
16
17 alpha=0.7
18 error_tol=1.e-5
19
20 # generate mesh
21
22 # print "\nGenerate mesh: finley.Rectangle(9,12,1)=>"
23 # mydomain=finley.Rectangle(140,140)
24
25 # print "\nGenerate mesh: finley.Rectangle(4,4,1)=>"
26 mydomain=finley.Rectangle(150,10,1)
27 # mydomain=finley.Rectangle(190,190,1)
28
29 print "\nGenerate mesh: finley.Rectangle(151,151,1)=>"
30 # mydomain=finley.Rectangle(151,151,1)
31 # mydomain=finley.Rectangle(128,128,1)
32
33 print "\nSetup domain and functions"
34 print "--------------------------"
35
36 print "e=Function(mydomain):"
37 e=Function(mydomain)
38
39 print "n=ContinuousFunction(mydomain):"
40 n=ContinuousFunction(mydomain)
41
42 # get handles to nodes and elements 1
43
44 print "\nGet handles to nodes and elements(1)=>"
45 print "--------------------------------------"
46
47 print "u_ex=Scalar(1,n,True):"
48 u_ex=Scalar(1,n,True)
49
50 print "x=e.getX():"
51 x=e.getX()
52
53 print "norm_u_ex=Lsup(u_ex):"
54 norm_u_ex=Lsup(u_ex)
55
56 print "\nGenerate a test solution (1)"
57 print "----------------------------"
58
59 print "mypde=LinearPDE( A=[[1.,0.8],[0.4,1.]], D=alpha, Y=alpha, domain=mydomain)"
60 mypde=LinearPDE(mydomain)
61 mypde.setDebugOn()
62 mypde.setValue(A=[[1.,0.1],[0.04,1.]],D=alpha,Y=alpha)
63
64 print "mypde.checkSymmetry()"
65 print mypde.checkSymmetry()
66
67 print "\nIterative Solver (1)=>"
68 mypde.setSolverMethod(mypde.BICGSTAB)
69 u_i=mypde.getSolution(verbose=True,iter_max=3000,preconditioner=mypde.ILU0)
70
71 print "\nDirect Solver (1)=>"
72 mypde.setSolverMethod(mypde.DIRECT)
73 u_d=mypde.getSolution(verbose=True)
74
75 print "\n***************************************************************"
76 error=u_ex-u_d
77 error_norm=Lsup(error)/norm_u_ex
78 print "norm of the error for direct solver is : ",error_norm
79 if error_norm > error_tol:
80 print "### error norm exceeded maximum tolerance ###"
81 sys.exit(1)
82 error=u_ex-u_i
83 error_norm=Lsup(error)/norm_u_ex
84 print "norm of the error for iterative solver is: ",error_norm
85 if error_norm > error_tol:
86 print "### error norm exceeded maximum tolerance ###"
87 sys.exit(1)
88 print "***************************************************************"
89
90 # get handles to nodes and elements 2
91
92 print "\nGet handles to nodes and elements(2)=>"
93 print "--------------------------------------"
94
95 print "x=n.getX():"
96 x=n.getX()
97
98 print " msk=whereZero(x[0])+whereZero(x[0]-1.)"
99 msk=whereZero(x[0])+whereZero(x[0]-1.)
100
101 print "mypde=LinearPDE(A=[[1.,0.],[0.,1.]],q=msk,r=u_ex)"
102 mypde=LinearPDE(mydomain)
103 mypde.setDebugOn()
104 mypde.setValue(A=[[1.,0.],[0.,1.]],q=msk,r=u_ex)
105
106 print "mypde.checkSymmetry()"
107 print mypde.checkSymmetry()
108
109 # generate a test solution 2
110
111 print "\nGenerate a test solution (2)"
112 print "----------------------------"
113
114 print "\nDirect Solver (2)=>"
115
116 mypde.setSymmetryOn()
117 mypde.setTolerance(1.e-13)
118
119 # mypde.setSymmetryOn() : is not woking yet!
120 mypde.setSolverMethod(mypde.DIRECT)
121 u_d=mypde.getSolution(verbose=True)
122
123 print "\nIterative Solver (2)=>"
124
125 mypde.setSolverMethod(mypde.ITERATIVE)
126 u_i=mypde.getSolution(verbose=True,iter_max=3000)
127
128 print "\n******************************************************************"
129 error=u_ex-u_d
130 error_norm=Lsup(error)/norm_u_ex
131 print "norm of the error for direct solver is : ",error_norm
132 if error_norm > error_tol:
133 print "### error norm exceeded maximum tolerance ###"
134 sys.exit(1)
135 error=u_ex-u_i
136 error_norm=Lsup(error)/norm_u_ex
137 print "norm of the error for iterative solver is: ",error_norm
138 if error_norm > error_tol:
139 print "### error norm exceeded maximum tolerance ###"
140 sys.exit(1)
141 print "******************************************************************"
142
143 print "\n-----"
144 print "Done."
145 print "-----"
146
147 stoptime = time.clock()
148 elapsed = stoptime - starttime
149 print "\nElapsed time: ", elapsed, "\n"
150
151 sys.exit(0)

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26