/[escript]/trunk-mpi-branch/finley/test/python/SimpleSolve.py
ViewVC logotype

Annotation of /trunk-mpi-branch/finley/test/python/SimpleSolve.py

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1016 - (hide annotations)
Thu Mar 8 06:31:28 2007 UTC (14 years, 6 months ago) by gross
File MIME type: text/x-python
File size: 4661 byte(s)
MPI version compiles and starts to run now. 
Important:

   * the mpi library needs to be shared.
   * the path needs to be added to LD_LIBRARY path.

The program stucks in the matrix assemblage.


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