test/python/SimpleSolve.py

# $Id$

import sys
import os
import unittest

from esys.escript import *
from esys.linearPDEs import *
from esys import finley

print "\nSimpleSolve.py"
print "--------------"

alpha=0.025

# generate mesh

# print "\nGenerate mesh: finley.Rectangle(9,12,1)=>"
# mydomain=finley.Rectangle(140,140)

# print "\nGenerate mesh: finley.Rectangle(4,4,1)=>"
# mydomain=finley.Rectangle(4,4,1)

print "\nGenerate mesh: finley.Rectangle(151,151,1)=>"
mydomain=finley.Rectangle(151,151,1)
# mydomain=finley.Rectangle(128,128,1)

print "\nSetup domain and functions"
print "--------------------------"

print "e=Function(mydomain):"
e=Function(mydomain)

print "n=ContinuousFunction(mydomain):"
n=ContinuousFunction(mydomain)

# get handles to nodes and elements 1

print "\nGet handles to nodes and elements(1)=>"
print "--------------------------------------"

print "u_ex=Scalar(1,n,True):"
u_ex=Scalar(1,n,True)

print "x=e.getX():"
x=e.getX()

print "norm_u_ex=u_ex.Lsup():"
norm_u_ex=u_ex.Lsup()

print "\nGenerate a test solution (1)"
print "----------------------------"

print "mypde=LinearPDE( A=[[1.,0.8],[0.4,1.]], D=alpha, Y=alpha, domain=mydomain)"
mypde=LinearPDE(mydomain)
mypde.setDebugOn()
mypde.setValue(A=[[1.,0.8],[0.4,1.]],D=alpha,Y=alpha)

print "mypde.checkSymmetry()"
print mypde.checkSymmetry()

print "\nIterative Solver (1)=>"
u_i=mypde.getSolution(preconditioner=ILU0,iter_max=3000)
# u_i=mypde.getSolution(iter_max=3000)


print "\nDirect Solver (1)=>"
mypde.setSolverMethod(DIRECT)
u_d=mypde.getSolution()


print "\n***************************************************************"
error=u_ex-u_d
print "norm of the error for direct solver is   : ",error.Lsup()/norm_u_ex
error=u_ex-u_i
print "norm of the error for iterative solver is: ",error.Lsup()/norm_u_ex
print "***************************************************************"

# get handles to nodes and elements 2

print "\nGet handles to nodes and elements(2)=>"
print "--------------------------------------"

print "x=n.getX():"
x=n.getX()

print "msk=x[0].whereZero()+(x[0]-1.).whereZero()"
msk=x[0].whereZero()+(x[0]-1.).whereZero()

print "mypde=LinearPDE(A=[[1.,0.],[0.,1.]],q=msk,r=u_ex)"
mypde=LinearPDE(mydomain)
mypde.setDebugOn()
mypde.setValue(A=[[1.,0.],[0.,1.]],q=msk,r=u_ex)

print "mypde.checkSymmetry()"
print mypde.checkSymmetry()

# generate a test solution 2

print "\nGenerate a test solution (2)"
print "----------------------------"

print "\nDirect Solver (2)=>"


# mypde.setSymmetryOn() : is not woking yet!
mypde.setSolverMethod(mypde.DIRECT)
u_d=mypde.getSolution()

print "\nIterative Solver (2)=>"

mypde.setSymmetryOn() 
mypde.setSolverMethod(DEFAULT_METHOD)
u_i=mypde.getSolution(iter_max=3000)

print "\n******************************************************************"
error=u_ex-u_d
print "norm of the error for direct solver is   : ",error.Lsup()/norm_u_ex
error=u_ex-u_i
print "norm of the error for iterative solver is: ",error.Lsup()/norm_u_ex
print "******************************************************************"

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