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