/[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 806 - (hide annotations)
Thu Aug 10 11:58:52 2006 UTC (14 years, 1 month ago) by gross
File MIME type: text/x-python
File size: 4600 byte(s)
Interface to the direct solver library UMLPACK is no implemented.


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 806 # mydomain=finley.Rectangle(5,5,1)
32 gross 777 mydomain=finley.Rectangle(500,500,1)
33 gross 584 # 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 gross 778 # set the direct solver switch
40     DIRECT=LinearPDE.DIRECT
41 gross 806 # DIRECT=LinearPDE.ITERATIVE # this will switch of the DIRECT solver call to avoid external library calls which may not be available everywhere
42 gross 778
43 jgs 82 print "\nSetup domain and functions"
44     print "--------------------------"
45    
46     print "e=Function(mydomain):"
47     e=Function(mydomain)
48    
49     print "n=ContinuousFunction(mydomain):"
50     n=ContinuousFunction(mydomain)
51    
52     # get handles to nodes and elements 1
53    
54     print "\nGet handles to nodes and elements(1)=>"
55     print "--------------------------------------"
56    
57     print "u_ex=Scalar(1,n,True):"
58     u_ex=Scalar(1,n,True)
59    
60     print "x=e.getX():"
61     x=e.getX()
62 gross 783 print "is x protected: ",x.isProtected()
63 jgs 82
64 gross 783
65 gross 298 print "norm_u_ex=Lsup(u_ex):"
66     norm_u_ex=Lsup(u_ex)
67 gross 783 print "is u_ex protected: ",u_ex.isProtected()
68 jgs 82
69 jgs 108 print "\nGenerate a test solution (1)"
70     print "----------------------------"
71    
72 jgs 102 print "mypde=LinearPDE( A=[[1.,0.8],[0.4,1.]], D=alpha, Y=alpha, domain=mydomain)"
73 jgs 104 mypde=LinearPDE(mydomain)
74 jgs 102 mypde.setDebugOn()
75 gross 432 mypde.setValue(A=[[1.,-0.001],[-0.001,1.]],D=alpha,Y=alpha)
76 jgs 82
77 jgs 108 print "mypde.checkSymmetry()"
78     print mypde.checkSymmetry()
79 jgs 82
80 jgs 102 print "\nIterative Solver (1)=>"
81 gross 453 # mypde.setSolverMethod(mypde.PRES20,preconditioner=mypde.ILU0)
82     mypde.setSolverMethod(mypde.BICGSTAB,preconditioner=mypde.JACOBI)
83 gross 431 u_i=mypde.getSolution(verbose=True,iter_max=3000)
84 jgs 102
85 jgs 100 print "\nDirect Solver (1)=>"
86 gross 778 mypde.setSolverMethod(DIRECT)
87 jgs 149 u_d=mypde.getSolution(verbose=True)
88 jgs 97
89 jgs 82 print "\n***************************************************************"
90     error=u_ex-u_d
91 gross 298 error_norm=Lsup(error)/norm_u_ex
92 jgs 147 print "norm of the error for direct solver is : ",error_norm
93     if error_norm > error_tol:
94     print "### error norm exceeded maximum tolerance ###"
95     sys.exit(1)
96 jgs 82 error=u_ex-u_i
97 gross 298 error_norm=Lsup(error)/norm_u_ex
98 jgs 147 print "norm of the error for iterative solver is: ",error_norm
99     if error_norm > error_tol:
100     print "### error norm exceeded maximum tolerance ###"
101     sys.exit(1)
102 jgs 82 print "***************************************************************"
103 gross 425 del mypde
104     print "***************************************************************"
105 jgs 82
106 gross 425
107 jgs 82 # get handles to nodes and elements 2
108    
109     print "\nGet handles to nodes and elements(2)=>"
110     print "--------------------------------------"
111    
112     print "x=n.getX():"
113     x=n.getX()
114    
115 gross 298 print " msk=whereZero(x[0])+whereZero(x[0]-1.)"
116     msk=whereZero(x[0])+whereZero(x[0]-1.)
117 jgs 82
118 jgs 102 print "mypde=LinearPDE(A=[[1.,0.],[0.,1.]],q=msk,r=u_ex)"
119 jgs 104 mypde=LinearPDE(mydomain)
120 jgs 108 mypde.setDebugOn()
121 gross 432 mypde.setValue(A=[[1.,0.0],[0.0,1.]],q=msk,r=u_ex)
122 jgs 82
123 jgs 108 print "mypde.checkSymmetry()"
124     print mypde.checkSymmetry()
125    
126 jgs 82 # generate a test solution 2
127    
128     print "\nGenerate a test solution (2)"
129     print "----------------------------"
130    
131     print "\nDirect Solver (2)=>"
132    
133 gross 404 mypde.setSymmetryOn()
134 jgs 123 mypde.setTolerance(1.e-13)
135 jgs 115
136 jgs 102 # mypde.setSymmetryOn() : is not woking yet!
137 gross 778 mypde.setSolverMethod(DIRECT)
138 jgs 149 u_d=mypde.getSolution(verbose=True)
139 jgs 82
140     print "\nIterative Solver (2)=>"
141    
142 gross 778 mypde.setSolverMethod(mypde.PCG)
143 jgs 149 u_i=mypde.getSolution(verbose=True,iter_max=3000)
144 jgs 82
145     print "\n******************************************************************"
146     error=u_ex-u_d
147 gross 298 error_norm=Lsup(error)/norm_u_ex
148 jgs 147 print "norm of the error for direct solver is : ",error_norm
149     if error_norm > error_tol:
150     print "### error norm exceeded maximum tolerance ###"
151     sys.exit(1)
152 jgs 82 error=u_ex-u_i
153 gross 298 error_norm=Lsup(error)/norm_u_ex
154 jgs 147 print "norm of the error for iterative solver is: ",error_norm
155 jgs 148 if error_norm > error_tol:
156 jgs 147 print "### error norm exceeded maximum tolerance ###"
157     sys.exit(1)
158 jgs 82 print "******************************************************************"
159    
160     print "\n-----"
161     print "Done."
162     print "-----"
163 jgs 147
164 jgs 153 stoptime = time.clock()
165     elapsed = stoptime - starttime
166     print "\nElapsed time: ", elapsed, "\n"
167    
168 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