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