/[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 153 - (hide annotations)
Tue Oct 25 01:51:20 2005 UTC (13 years, 11 months ago) by jgs
Original Path: trunk/esys2/finley/test/python/SimpleSolve.py
File MIME type: text/x-python
File size: 3886 byte(s)
Merge of development branch dev-02 back to main trunk on 2005-10-25

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 jgs 150 # mydomain=finley.Rectangle(190,190,1)
28 jgs 82
29     print "\nGenerate mesh: finley.Rectangle(151,151,1)=>"
30 jgs 148 # mydomain=finley.Rectangle(151,151,1)
31 jgs 115 # mydomain=finley.Rectangle(128,128,1)
32 jgs 82
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=u_ex.Lsup():"
54     norm_u_ex=u_ex.Lsup()
55    
56 jgs 108 print "\nGenerate a test solution (1)"
57     print "----------------------------"
58    
59 jgs 102 print "mypde=LinearPDE( A=[[1.,0.8],[0.4,1.]], D=alpha, Y=alpha, domain=mydomain)"
60 jgs 104 mypde=LinearPDE(mydomain)
61 jgs 102 mypde.setDebugOn()
62 jgs 150 mypde.setValue(A=[[1.,0.1],[0.04,1.]],D=alpha,Y=alpha)
63 jgs 82
64 jgs 108 print "mypde.checkSymmetry()"
65     print mypde.checkSymmetry()
66 jgs 82
67 jgs 102 print "\nIterative Solver (1)=>"
68 jgs 150 mypde.setSolverMethod(mypde.BICGSTAB)
69 jgs 149 u_i=mypde.getSolution(verbose=True,iter_max=3000,preconditioner=mypde.ILU0)
70 jgs 102
71 jgs 100 print "\nDirect Solver (1)=>"
72 jgs 148 mypde.setSolverMethod(mypde.DIRECT)
73 jgs 149 u_d=mypde.getSolution(verbose=True)
74 jgs 97
75 jgs 82 print "\n***************************************************************"
76     error=u_ex-u_d
77 jgs 147 error_norm=error.Lsup()/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 jgs 82 error=u_ex-u_i
83 jgs 147 error_norm=error.Lsup()/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 jgs 82 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=x[0].whereZero()+(x[0]-1.).whereZero()"
99     msk=x[0].whereZero()+(x[0]-1.).whereZero()
100    
101 jgs 102 print "mypde=LinearPDE(A=[[1.,0.],[0.,1.]],q=msk,r=u_ex)"
102 jgs 104 mypde=LinearPDE(mydomain)
103 jgs 108 mypde.setDebugOn()
104 jgs 104 mypde.setValue(A=[[1.,0.],[0.,1.]],q=msk,r=u_ex)
105 jgs 82
106 jgs 108 print "mypde.checkSymmetry()"
107     print mypde.checkSymmetry()
108    
109 jgs 82 # generate a test solution 2
110    
111     print "\nGenerate a test solution (2)"
112     print "----------------------------"
113    
114     print "\nDirect Solver (2)=>"
115    
116 jgs 123 #mypde.setSymmetryOn()
117     mypde.setTolerance(1.e-13)
118 jgs 115
119 jgs 102 # mypde.setSymmetryOn() : is not woking yet!
120 jgs 104 mypde.setSolverMethod(mypde.DIRECT)
121 jgs 149 u_d=mypde.getSolution(verbose=True)
122 jgs 82
123     print "\nIterative Solver (2)=>"
124    
125 jgs 108 mypde.setSymmetryOn()
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 jgs 147 error_norm=error.Lsup()/norm_u_ex
132     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 jgs 147 error_norm=error.Lsup()/norm_u_ex
138     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