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

1 jgs 82 # $Id$
2    
3 elspeth 617 __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 jgs 82 import sys
9     import os
10     import unittest
11 jgs 153 import time
12 jgs 82
13 jgs 149 from esys.escript import *
14     from esys.escript.linearPDEs import *
15     from esys import finley
16 jgs 82
17 jgs 153 starttime = time.clock()
18    
19 jgs 82 print "\nSimpleSolve.py"
20     print "--------------"
21    
22 jgs 150 alpha=0.7
23 jgs 148 error_tol=1.e-5
24 jgs 82
25     # generate mesh
26    
27 jgs 102 # print "\nGenerate mesh: finley.Rectangle(9,12,1)=>"
28 jgs 115 # mydomain=finley.Rectangle(140,140)
29 jgs 82
30 jgs 102 # print "\nGenerate mesh: finley.Rectangle(4,4,1)=>"
31 gross 584 mydomain=finley.Rectangle(50,50,1)
32     #mydomain=finley.Rectangle(500,500,1)
33     # mydomain=finley.Rectangle(150,150,1)
34 jgs 82
35     print "\nGenerate mesh: finley.Rectangle(151,151,1)=>"
36 jgs 148 # mydomain=finley.Rectangle(151,151,1)
37 jgs 115 # mydomain=finley.Rectangle(128,128,1)
38 jgs 82
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 gross 298 print "norm_u_ex=Lsup(u_ex):"
60     norm_u_ex=Lsup(u_ex)
61 jgs 82
62 jgs 108 print "\nGenerate a test solution (1)"
63     print "----------------------------"
64    
65 jgs 102 print "mypde=LinearPDE( A=[[1.,0.8],[0.4,1.]], D=alpha, Y=alpha, domain=mydomain)"
66 jgs 104 mypde=LinearPDE(mydomain)
67 jgs 102 mypde.setDebugOn()
68 gross 432 mypde.setValue(A=[[1.,-0.001],[-0.001,1.]],D=alpha,Y=alpha)
69 jgs 82
70 jgs 108 print "mypde.checkSymmetry()"
71     print mypde.checkSymmetry()
72 jgs 82
73 jgs 102 print "\nIterative Solver (1)=>"
74 gross 453 # mypde.setSolverMethod(mypde.PRES20,preconditioner=mypde.ILU0)
75     mypde.setSolverMethod(mypde.BICGSTAB,preconditioner=mypde.JACOBI)
76 gross 431 u_i=mypde.getSolution(verbose=True,iter_max=3000)
77 jgs 102
78 jgs 100 print "\nDirect Solver (1)=>"
79 jgs 148 mypde.setSolverMethod(mypde.DIRECT)
80 jgs 149 u_d=mypde.getSolution(verbose=True)
81 jgs 97
82 jgs 82 print "\n***************************************************************"
83     error=u_ex-u_d
84 gross 298 error_norm=Lsup(error)/norm_u_ex
85 jgs 147 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 jgs 82 error=u_ex-u_i
90 gross 298 error_norm=Lsup(error)/norm_u_ex
91 jgs 147 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 jgs 82 print "***************************************************************"
96 gross 425 del mypde
97     print "***************************************************************"
98 jgs 82
99 gross 425
100 jgs 82 # 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 gross 298 print " msk=whereZero(x[0])+whereZero(x[0]-1.)"
109     msk=whereZero(x[0])+whereZero(x[0]-1.)
110 jgs 82
111 jgs 102 print "mypde=LinearPDE(A=[[1.,0.],[0.,1.]],q=msk,r=u_ex)"
112 jgs 104 mypde=LinearPDE(mydomain)
113 jgs 108 mypde.setDebugOn()
114 gross 432 mypde.setValue(A=[[1.,0.0],[0.0,1.]],q=msk,r=u_ex)
115 jgs 82
116 jgs 108 print "mypde.checkSymmetry()"
117     print mypde.checkSymmetry()
118    
119 jgs 82 # generate a test solution 2
120    
121     print "\nGenerate a test solution (2)"
122     print "----------------------------"
123    
124     print "\nDirect Solver (2)=>"
125    
126 gross 404 mypde.setSymmetryOn()
127 jgs 123 mypde.setTolerance(1.e-13)
128 jgs 115
129 jgs 102 # mypde.setSymmetryOn() : is not woking yet!
130 jgs 104 mypde.setSolverMethod(mypde.DIRECT)
131 jgs 149 u_d=mypde.getSolution(verbose=True)
132 jgs 82
133     print "\nIterative Solver (2)=>"
134    
135 jgs 150 mypde.setSolverMethod(mypde.ITERATIVE)
136 jgs 149 u_i=mypde.getSolution(verbose=True,iter_max=3000)
137 jgs 82
138     print "\n******************************************************************"
139     error=u_ex-u_d
140 gross 298 error_norm=Lsup(error)/norm_u_ex
141 jgs 147 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 jgs 82 error=u_ex-u_i
146 gross 298 error_norm=Lsup(error)/norm_u_ex
147 jgs 147 print "norm of the error for iterative solver is: ",error_norm
148 jgs 148 if error_norm > error_tol:
149 jgs 147 print "### error norm exceeded maximum tolerance ###"
150     sys.exit(1)
151 jgs 82 print "******************************************************************"
152    
153     print "\n-----"
154     print "Done."
155     print "-----"
156 jgs 147
157 jgs 153 stoptime = time.clock()
158     elapsed = stoptime - starttime
159     print "\nElapsed time: ", elapsed, "\n"
160    
161 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