/[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 150 - (hide annotations)
Thu Sep 15 03:44:45 2005 UTC (14 years, 6 months ago) by jgs
Original Path: trunk/esys2/finley/test/python/SimpleSolve.py
File MIME type: text/x-python
File size: 3751 byte(s)
Merge of development branch dev-02 back to main trunk on 2005-09-15

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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26