/[escript]/trunk/finley/test/python/blocktest.py
ViewVC logotype

Contents of /trunk/finley/test/python/blocktest.py

Parent Directory Parent Directory | Revision Log Revision Log


Revision 6651 - (show annotations)
Wed Feb 7 02:12:08 2018 UTC (20 months, 1 week ago) by jfenwick
File MIME type: text/x-python
File size: 3534 byte(s)
Make everyone sad by touching all the files

Copyright dates update

1 #
2 # this script is testing block solvers for PDEs
3 #
4 #
5 # - u_{j,ii} + b*u_j+ a*sum_{k<>j} (u_j-u_k) = F_j
6 #
7 # where a controls the degree of coupling and b the degree of diagonal dominance.
8 # a and b may have any value.
9 #
10 # The domain needs to be a unit square or cube with any type of mesh
11 #
12 #
13 ##############################################################################
14 #
15 # Copyright (c) 2003-2018 by The University of Queensland
16 # http://www.uq.edu.au
17 #
18 # Primary Business: Queensland, Australia
19 # Licensed under the Apache License, version 2.0
20 # http://www.apache.org/licenses/LICENSE-2.0
21 #
22 # Development until 2012 by Earth Systems Science Computational Center (ESSCC)
23 # Development 2012-2013 by School of Earth Sciences
24 # Development from 2014 by Centre for Geoscience Computing (GeoComp)
25 #
26 ##############################################################################
27
28 from __future__ import print_function, division
29
30 __copyright__="""Copyright (c) 2003-2018 by The University of Queensland
31 http://www.uq.edu.au
32 Primary Business: Queensland, Australia"""
33 __license__="""Licensed under the Apache License, version 2.0
34 http://www.apache.org/licenses/LICENSE-2.0"""
35 __url__="https://launchpad.net/escript-finley"
36
37 from esys.escript import *
38 from esys.escript.linearPDEs import LinearPDE
39 import esys.finley as finley
40
41 TOL=1.e-8
42
43 def runTest(dom, n=1, a=0, b=0):
44 print("================== TEST : n= %s a=%s b=%s ================"%(n,a,b))
45 DIM=dom.getDim()
46 normal=dom.getNormal()
47 mypde=LinearPDE(dom, numEquations=n, numSolutions=n)
48 x=dom.getX()
49 A=mypde.createCoefficient("A")
50 D=mypde.createCoefficient("D")
51 Y=mypde.createCoefficient("Y")
52 y=mypde.createCoefficient("y")
53 q=mypde.createCoefficient("q")
54 if n==1:
55 u_ref=Scalar(0.,Solution(dom))
56 for j in range(DIM):
57 q+=whereZero(x[j])
58 A[j,j]=1
59 y+=sin(sqrt(2.))*normal[0]
60 Y+=b*x[0]*sin(sqrt(2.))
61 D+=b
62 u_ref=x[0]*sin(sqrt(2.))
63 else:
64 u_ref=Data(0.,(n,),Solution(dom))
65 for i in range(n):
66 for j in range(DIM):
67 q[i]+=whereZero(x[j])
68 A[i,j,i,j]=1
69 if j == i%DIM: y[i]+=sin(i+sqrt(2.))*normal[j]
70 for j in range(n):
71 if i==j:
72 D[i,i]=b+(n-1)*a
73 Y[i]+=b*x[i%DIM]*sin(i+sqrt(2.))
74 else:
75 D[i,j]=-a
76 Y[i]+=a*(x[i%DIM]*sin(i+sqrt(2.))-x[j%DIM]*sin(j+sqrt(2.)))
77 u_ref[i]=x[i%DIM]*sin(i+sqrt(2.))
78
79 # - u_{j,ii} + b*u_i+ a*sum_{k<>j} (u_i-u_k) = F_j
80
81 mypde.setValue(A=A, D=D, y=y, q=q, r=u_ref, Y=Y)
82 mypde.getSolverOptions().setVerbosityOn()
83 mypde.getSolverOptions().setTolerance(TOL)
84 mypde.setSymmetryOn()
85 u=mypde.getSolution()
86 error=Lsup(u-u_ref)/Lsup(u_ref)
87 print("error = ",error)
88 if error > 10*TOL: print("XXXXXXXXXX Error to large ! XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX")
89
90 domain=finley.Rectangle(10,20,1,l0=1.,l1=1.0)
91 # or Brick or ReadMesh
92 runTest(dom=domain, n=1, b=0)
93 runTest(dom=domain, n=1, b=5)
94 runTest(dom=domain, n=1, b=50)
95
96 runTest(dom=domain, n=2, a=0, b=0)
97 runTest(dom=domain, n=2, a=5, b=0)
98 runTest(dom=domain, n=2, a=50, b=0)
99
100 runTest(dom=domain, n=2, a=0, b=10)
101 runTest(dom=domain, n=2, a=5, b=10)
102 runTest(dom=domain, n=2, a=50, b=10)
103
104 runTest(dom=domain, n=3, a=0, b=0)
105 runTest(dom=domain, n=3, a=5, b=0)
106 runTest(dom=domain, n=3, a=50, b=0)
107
108 runTest(dom=domain, n=3, a=0, b=10)
109 runTest(dom=domain, n=3, a=5, b=10)
110 runTest(dom=domain, n=3, a=50, b=10)

Properties

Name Value
svn:executable *

  ViewVC Help
Powered by ViewVC 1.1.26