/[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 432 - (hide annotations)
Fri Jan 13 07:38:54 2006 UTC (13 years, 10 months ago) by gross
File MIME type: text/x-python
File size: 3973 byte(s)
some fixes for openmp
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 gross 432 mydomain=finley.Rectangle(60,40,1)
27     # mydomain=finley.Rectangle(250,250,1)
28     mydomain=finley.Rectangle(100,100,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 gross 432 mypde.setValue(A=[[1.,-0.001],[-0.001,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 gross 432 mypde.setSolverMethod(mypde.PRES20,preconditioner=mypde.ILU0)
70 gross 431 u_i=mypde.getSolution(verbose=True,iter_max=3000)
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 gross 425 del mypde
91     print "***************************************************************"
92 jgs 82
93 gross 425
94 jgs 82 # get handles to nodes and elements 2
95    
96     print "\nGet handles to nodes and elements(2)=>"
97     print "--------------------------------------"
98    
99     print "x=n.getX():"
100     x=n.getX()
101    
102 gross 298 print " msk=whereZero(x[0])+whereZero(x[0]-1.)"
103     msk=whereZero(x[0])+whereZero(x[0]-1.)
104 jgs 82
105 jgs 102 print "mypde=LinearPDE(A=[[1.,0.],[0.,1.]],q=msk,r=u_ex)"
106 jgs 104 mypde=LinearPDE(mydomain)
107 jgs 108 mypde.setDebugOn()
108 gross 432 mypde.setValue(A=[[1.,0.0],[0.0,1.]],q=msk,r=u_ex)
109 jgs 82
110 jgs 108 print "mypde.checkSymmetry()"
111     print mypde.checkSymmetry()
112    
113 jgs 82 # generate a test solution 2
114    
115     print "\nGenerate a test solution (2)"
116     print "----------------------------"
117    
118     print "\nDirect Solver (2)=>"
119    
120 gross 404 mypde.setSymmetryOn()
121 jgs 123 mypde.setTolerance(1.e-13)
122 jgs 115
123 jgs 102 # mypde.setSymmetryOn() : is not woking yet!
124 jgs 104 mypde.setSolverMethod(mypde.DIRECT)
125 jgs 149 u_d=mypde.getSolution(verbose=True)
126 jgs 82
127     print "\nIterative Solver (2)=>"
128    
129 jgs 150 mypde.setSolverMethod(mypde.ITERATIVE)
130 jgs 149 u_i=mypde.getSolution(verbose=True,iter_max=3000)
131 jgs 82
132     print "\n******************************************************************"
133     error=u_ex-u_d
134 gross 298 error_norm=Lsup(error)/norm_u_ex
135 jgs 147 print "norm of the error for direct solver is : ",error_norm
136     if error_norm > error_tol:
137     print "### error norm exceeded maximum tolerance ###"
138     sys.exit(1)
139 jgs 82 error=u_ex-u_i
140 gross 298 error_norm=Lsup(error)/norm_u_ex
141 jgs 147 print "norm of the error for iterative solver is: ",error_norm
142 jgs 148 if error_norm > error_tol:
143 jgs 147 print "### error norm exceeded maximum tolerance ###"
144     sys.exit(1)
145 jgs 82 print "******************************************************************"
146    
147     print "\n-----"
148     print "Done."
149     print "-----"
150 jgs 147
151 jgs 153 stoptime = time.clock()
152     elapsed = stoptime - starttime
153     print "\nElapsed time: ", elapsed, "\n"
154    
155 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