1 |
# $Id$ |
2 |
|
3 |
import sys |
4 |
import os |
5 |
import unittest |
6 |
import time |
7 |
|
8 |
from esys.escript import * |
9 |
from esys.escript.linearPDEs import * |
10 |
from esys import finley |
11 |
|
12 |
starttime = time.clock() |
13 |
|
14 |
print "\nSimpleSolve.py" |
15 |
print "--------------" |
16 |
|
17 |
alpha=0.7 |
18 |
error_tol=1.e-5 |
19 |
|
20 |
# generate mesh |
21 |
|
22 |
# print "\nGenerate mesh: finley.Rectangle(9,12,1)=>" |
23 |
# mydomain=finley.Rectangle(140,140) |
24 |
|
25 |
# print "\nGenerate mesh: finley.Rectangle(4,4,1)=>" |
26 |
mydomain=finley.Rectangle(150,10,1) |
27 |
# mydomain=finley.Rectangle(190,190,1) |
28 |
|
29 |
print "\nGenerate mesh: finley.Rectangle(151,151,1)=>" |
30 |
# mydomain=finley.Rectangle(151,151,1) |
31 |
# mydomain=finley.Rectangle(128,128,1) |
32 |
|
33 |
print "\nSetup domain and functions" |
34 |
print "--------------------------" |
35 |
|
36 |
print "e=Function(mydomain):" |
37 |
e=Function(mydomain) |
38 |
|
39 |
print "n=ContinuousFunction(mydomain):" |
40 |
n=ContinuousFunction(mydomain) |
41 |
|
42 |
# get handles to nodes and elements 1 |
43 |
|
44 |
print "\nGet handles to nodes and elements(1)=>" |
45 |
print "--------------------------------------" |
46 |
|
47 |
print "u_ex=Scalar(1,n,True):" |
48 |
u_ex=Scalar(1,n,True) |
49 |
|
50 |
print "x=e.getX():" |
51 |
x=e.getX() |
52 |
|
53 |
print "norm_u_ex=Lsup(u_ex):" |
54 |
norm_u_ex=Lsup(u_ex) |
55 |
|
56 |
print "\nGenerate a test solution (1)" |
57 |
print "----------------------------" |
58 |
|
59 |
print "mypde=LinearPDE( A=[[1.,0.8],[0.4,1.]], D=alpha, Y=alpha, domain=mydomain)" |
60 |
mypde=LinearPDE(mydomain) |
61 |
mypde.setDebugOn() |
62 |
mypde.setValue(A=[[1.,0.1],[0.04,1.]],D=alpha,Y=alpha) |
63 |
|
64 |
print "mypde.checkSymmetry()" |
65 |
print mypde.checkSymmetry() |
66 |
|
67 |
print "\nIterative Solver (1)=>" |
68 |
mypde.setSolverMethod(mypde.BICGSTAB) |
69 |
u_i=mypde.getSolution(verbose=True,iter_max=3000,preconditioner=mypde.ILU0) |
70 |
|
71 |
print "\nDirect Solver (1)=>" |
72 |
mypde.setSolverMethod(mypde.DIRECT) |
73 |
u_d=mypde.getSolution(verbose=True) |
74 |
|
75 |
print "\n***************************************************************" |
76 |
error=u_ex-u_d |
77 |
error_norm=Lsup(error)/norm_u_ex |
78 |
print "norm of the error for direct solver is : ",error_norm |
79 |
if error_norm > error_tol: |
80 |
print "### error norm exceeded maximum tolerance ###" |
81 |
sys.exit(1) |
82 |
error=u_ex-u_i |
83 |
error_norm=Lsup(error)/norm_u_ex |
84 |
print "norm of the error for iterative solver is: ",error_norm |
85 |
if error_norm > error_tol: |
86 |
print "### error norm exceeded maximum tolerance ###" |
87 |
sys.exit(1) |
88 |
print "***************************************************************" |
89 |
|
90 |
# get handles to nodes and elements 2 |
91 |
|
92 |
print "\nGet handles to nodes and elements(2)=>" |
93 |
print "--------------------------------------" |
94 |
|
95 |
print "x=n.getX():" |
96 |
x=n.getX() |
97 |
|
98 |
print " msk=whereZero(x[0])+whereZero(x[0]-1.)" |
99 |
msk=whereZero(x[0])+whereZero(x[0]-1.) |
100 |
|
101 |
print "mypde=LinearPDE(A=[[1.,0.],[0.,1.]],q=msk,r=u_ex)" |
102 |
mypde=LinearPDE(mydomain) |
103 |
mypde.setDebugOn() |
104 |
mypde.setValue(A=[[1.,0.],[0.,1.]],q=msk,r=u_ex) |
105 |
|
106 |
print "mypde.checkSymmetry()" |
107 |
print mypde.checkSymmetry() |
108 |
|
109 |
# generate a test solution 2 |
110 |
|
111 |
print "\nGenerate a test solution (2)" |
112 |
print "----------------------------" |
113 |
|
114 |
print "\nDirect Solver (2)=>" |
115 |
|
116 |
mypde.setSymmetryOn() |
117 |
mypde.setTolerance(1.e-13) |
118 |
|
119 |
# mypde.setSymmetryOn() : is not woking yet! |
120 |
mypde.setSolverMethod(mypde.DIRECT) |
121 |
u_d=mypde.getSolution(verbose=True) |
122 |
|
123 |
print "\nIterative Solver (2)=>" |
124 |
|
125 |
mypde.setSolverMethod(mypde.ITERATIVE) |
126 |
u_i=mypde.getSolution(verbose=True,iter_max=3000) |
127 |
|
128 |
print "\n******************************************************************" |
129 |
error=u_ex-u_d |
130 |
error_norm=Lsup(error)/norm_u_ex |
131 |
print "norm of the error for direct solver is : ",error_norm |
132 |
if error_norm > error_tol: |
133 |
print "### error norm exceeded maximum tolerance ###" |
134 |
sys.exit(1) |
135 |
error=u_ex-u_i |
136 |
error_norm=Lsup(error)/norm_u_ex |
137 |
print "norm of the error for iterative solver is: ",error_norm |
138 |
if error_norm > error_tol: |
139 |
print "### error norm exceeded maximum tolerance ###" |
140 |
sys.exit(1) |
141 |
print "******************************************************************" |
142 |
|
143 |
print "\n-----" |
144 |
print "Done." |
145 |
print "-----" |
146 |
|
147 |
stoptime = time.clock() |
148 |
elapsed = stoptime - starttime |
149 |
print "\nElapsed time: ", elapsed, "\n" |
150 |
|
151 |
sys.exit(0) |