test/python/SimpleSolve.py

# $Id$

import sys
import os
import unittest

esys_root=os.getenv('ESYS_ROOT')
sys.path.append(esys_root+'/finley/lib')
sys.path.append(esys_root+'/escript/lib')
sys.path.append(esys_root+'/escript/py_src')

from escript import *
from util import *
from linearPDE import *
import finley

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

my_options = {
          "verbose" : True,
          "reordering" : NO_REORDERING,
          "tolerance" : 1.E-8,
          "final_residual" : 0.,
          "iterative_method" : PCG,
          "preconditioner" : JACOBI,
          "iter_max" :  5000,
          "iter" : 0,
          "drop_tolerance" : 1.10,
          "drop_storage" : 2.
}

alpha=0.01

print "\nOptions: ", my_options

# generate mesh

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

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)

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 "mypde=linearPDE( A=[[1.,0.7],[0.7,1.]], D=alpha, Y=alpha, domain=mydomain)"
mypde=linearPDE(A=[[1.,0.7],[0.7,1.]],D=alpha,Y=alpha,domain=mydomain)
#mypde=linearPDE(A=[[1.,0.],[0.,1.]],D=alpha,Y=alpha,domain=mydomain)

# generate a test solution 1

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

print "\nDirect Solver (1)=>"

u_d=mypde.getSolution(iterative=False,**my_options)

print "\nIterative Solver (1)=>"

u_i=mypde.getSolution(iterative=True,**my_options)

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(A=[[1.,0.],[0.,1.]],q=msk,r=u_ex)

# generate a test solution 2

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

print "\nDirect Solver (2)=>"

u_d=mypde.getSolution(iterative=False,**my_options)

print "\nIterative Solver (2)=>"

u_i=mypde.getSolution(iterative=True,**my_options)

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	jgs	82	# $Id$
2
3			import sys
4			import os
5			import unittest
6
7			esys_root=os.getenv('ESYS_ROOT')
8			sys.path.append(esys_root+'/finley/lib')
9			sys.path.append(esys_root+'/escript/lib')
10			sys.path.append(esys_root+'/escript/py_src')
11
12			from escript import *
13			from util import *
14			from linearPDE import *
15			import finley
16
17			print "\nSimpleSolve.py"
18			print "--------------"
19
20			my_options = {
21			"verbose" : True,
22			"reordering" : NO_REORDERING,
23			"tolerance" : 1.E-8,
24			"final_residual" : 0.,
25			"iterative_method" : PCG,
26			"preconditioner" : JACOBI,
27			"iter_max" : 5000,
28			"iter" : 0,
29			"drop_tolerance" : 1.10,
30			"drop_storage" : 2.
31			}
32
33			alpha=0.01
34
35			print "\nOptions: ", my_options
36
37			# generate mesh
38
39			print "\nGenerate mesh: finley.Rectangle(9,12,1)=>"
40			mydomain=finley.Rectangle(9,12,1)
41
42			print "\nGenerate mesh: finley.Rectangle(4,4,1)=>"
43			mydomain=finley.Rectangle(4,4,1)
44
45			print "\nGenerate mesh: finley.Rectangle(151,151,1)=>"
46			mydomain=finley.Rectangle(151,151,1)
47
48			print "\nSetup domain and functions"
49			print "--------------------------"
50
51			print "e=Function(mydomain):"
52			e=Function(mydomain)
53
54			print "n=ContinuousFunction(mydomain):"
55			n=ContinuousFunction(mydomain)
56
57			# get handles to nodes and elements 1
58
59			print "\nGet handles to nodes and elements(1)=>"
60			print "--------------------------------------"
61
62			print "u_ex=Scalar(1,n,True):"
63			u_ex=Scalar(1,n,True)
64
65			print "x=e.getX():"
66			x=e.getX()
67
68			print "norm_u_ex=u_ex.Lsup():"
69			norm_u_ex=u_ex.Lsup()
70
71			print "mypde=linearPDE( A=[[1.,0.7],[0.7,1.]], D=alpha, Y=alpha, domain=mydomain)"
72			mypde=linearPDE(A=[[1.,0.7],[0.7,1.]],D=alpha,Y=alpha,domain=mydomain)
73			#mypde=linearPDE(A=[[1.,0.],[0.,1.]],D=alpha,Y=alpha,domain=mydomain)
74
75			# generate a test solution 1
76
77			print "\nGenerate a test solution (1)"
78			print "----------------------------"
79
80			print "\nDirect Solver (1)=>"
81
82			u_d=mypde.getSolution(iterative=False,**my_options)
83
84			print "\nIterative Solver (1)=>"
85
86			u_i=mypde.getSolution(iterative=True,**my_options)
87
88			print "\n***************************************************************"
89			error=u_ex-u_d
90			print "norm of the error for direct solver is : ",error.Lsup()/norm_u_ex
91			error=u_ex-u_i
92			print "norm of the error for iterative solver is: ",error.Lsup()/norm_u_ex
93			print "***************************************************************"
94
95			# get handles to nodes and elements 2
96
97			print "\nGet handles to nodes and elements(2)=>"
98			print "--------------------------------------"
99
100			print "x=n.getX():"
101			x=n.getX()
102
103			print "msk=x[0].whereZero()+(x[0]-1.).whereZero()"
104			msk=x[0].whereZero()+(x[0]-1.).whereZero()
105
106			print "mypde=linearPDE(A=[[1.,0.],[0.,1.]],q=msk,r=u_ex)"
107			mypde=linearPDE(A=[[1.,0.],[0.,1.]],q=msk,r=u_ex)
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			u_d=mypde.getSolution(iterative=False,**my_options)
117
118			print "\nIterative Solver (2)=>"
119
120			u_i=mypde.getSolution(iterative=True,**my_options)
121
122			print "\n******************************************************************"
123			error=u_ex-u_d
124			print "norm of the error for direct solver is : ",error.Lsup()/norm_u_ex
125			error=u_ex-u_i
126			print "norm of the error for iterative solver is: ",error.Lsup()/norm_u_ex
127			print "******************************************************************"
128
129			print "\n-----"
130			print "Done."
131			print "-----"
Name	Value
svn:eol-style	native
svn:keywords	Author Date Id Revision