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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1016 - (show annotations)
Thu Mar 8 06:31:28 2007 UTC (14 years, 4 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 # $Id$
2
3 __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 import sys
9 import os
10 import unittest
11 import time
12
13 from esys.escript import *
14 from esys.escript.linearPDEs import *
15 from esys import finley
16
17 starttime = time.time()
18
19 print "\nSimpleSolve.py"
20 print "--------------"
21
22 alpha=0.7
23 error_tol=1.e-5
24
25 # generate mesh
26
27 # print "\nGenerate mesh: finley.Rectangle(9,12,1)=>"
28 # mydomain=finley.Rectangle(140,140)
29
30 # print "\nGenerate mesh: finley.Rectangle(4,4,1)=>"
31 # mydomain=finley.Rectangle(5,5,1)
32 print "A"
33 mydomain=finley.Rectangle(500,500,1)
34 print "B"
35 mydomain=finley.Rectangle(150,150,1)
36 print "C"
37
38 print "\nGenerate mesh: finley.Rectangle(151,151,1)=>"
39 # mydomain=finley.Rectangle(151,151,1)
40 # mydomain=finley.Rectangle(128,128,1)
41
42 # set the direct solver switch
43 DIRECT=LinearPDE.DIRECT
44 # DIRECT=LinearPDE.ITERATIVE # this will switch of the DIRECT solver call to avoid external library calls which may not be available everywhere
45
46 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 print "is x protected: ",x.isProtected()
66
67
68 print "norm_u_ex=Lsup(u_ex):"
69 norm_u_ex=Lsup(u_ex)
70 print "is u_ex protected: ",u_ex.isProtected()
71
72 print "\nGenerate a test solution (1)"
73 print "----------------------------"
74
75 print "mypde=LinearPDE( A=[[1.,0.8],[0.4,1.]], D=alpha, Y=alpha, domain=mydomain)"
76 mypde=LinearPDE(mydomain)
77 mypde.setDebugOn()
78 mypde.setValue(A=[[1.,-0.001],[-0.001,1.]],D=alpha,Y=alpha)
79
80 print "mypde.checkSymmetry()"
81 print mypde.checkSymmetry()
82
83 print "\nIterative Solver (1)=>"
84 # mypde.setSolverMethod(mypde.PRES20,preconditioner=mypde.ILU0)
85 mypde.setSolverMethod(mypde.BICGSTAB,preconditioner=mypde.JACOBI)
86 u_i=mypde.getSolution(verbose=True,iter_max=3000)
87
88 print "\nDirect Solver (1)=>"
89 mypde.setSolverMethod(DIRECT)
90 u_d=mypde.getSolution(verbose=True)
91
92 print "\n***************************************************************"
93 error=u_ex-u_d
94 error_norm=Lsup(error)/norm_u_ex
95 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 error=u_ex-u_i
100 error_norm=Lsup(error)/norm_u_ex
101 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 print "***************************************************************"
106 # del mypde
107 del u_i
108 # releaseUnusedMemory()
109 print "***************************************************************"
110
111
112
113 # 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 print " msk=whereZero(x[0])+whereZero(x[0]-1.)"
122 msk=whereZero(x[0])+whereZero(x[0]-1.)
123
124 print "mypde=LinearPDE(A=[[1.,0.],[0.,1.]],q=msk,r=u_ex)"
125 mypde=LinearPDE(mydomain)
126 mypde.setDebugOn()
127 mypde.setValue(A=[[1.,0.0],[0.0,1.]],q=msk,r=u_ex)
128
129 print "mypde.checkSymmetry()"
130 print mypde.checkSymmetry()
131
132 # generate a test solution 2
133
134 print "\nGenerate a test solution (2)"
135 print "----------------------------"
136
137 print "\nDirect Solver (2)=>"
138
139 mypde.setSymmetryOn()
140 mypde.setTolerance(1.e-13)
141
142 # mypde.setSymmetryOn() : is not woking yet!
143 mypde.setSolverMethod(DIRECT)
144 u_d=mypde.getSolution(verbose=True)
145
146 print "\nIterative Solver (2)=>"
147
148 mypde.setSolverMethod(mypde.PCG)
149 u_i=mypde.getSolution(verbose=True,iter_max=3000)
150
151 print "\n******************************************************************"
152 error=u_ex-u_d
153 error_norm=Lsup(error)/norm_u_ex
154 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 error=u_ex-u_i
159 error_norm=Lsup(error)/norm_u_ex
160 print "norm of the error for iterative solver is: ",error_norm
161 if error_norm > error_tol:
162 print "### error norm exceeded maximum tolerance ###"
163 sys.exit(1)
164 print "******************************************************************"
165
166 print "\n-----"
167 print "Done."
168 print "-----"
169
170 stoptime = time.time()
171 elapsed = stoptime - starttime
172 print "\nElapsed time: ", elapsed, "\n"
173
174 sys.exit(0)

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26