# Contents of /trunk/finley/test/python/SolveTest.py

Revision 102 - (show annotations)
Wed Dec 15 07:08:39 2004 UTC (15 years, 3 months ago) by jgs
Original Path: trunk/esys2/finley/test/python/SolveTest.py
File MIME type: text/x-python
File size: 6002 byte(s)
```*** empty log message ***

```
 1 # \$Id\$ 2 3 """General test environment to test the solvers for scalar and vector equations 4 5 test parameters are 6 7 numDim = spatial dimension 8 totalNumElem = number of func in each direction 9 problem = solveScalar,solveVector 10 11 solver_method = true/false 12 len_x0 = length of the domain in x0 direction (number of func in x0 is round(totalNumElem*len_x0) ) 13 alpha = a parameter of the PDE (not well defined yet) 14 15 """ 16 17 import sys 18 import os 19 20 esys_root=os.getenv('ESYS_ROOT') 21 sys.path.append(esys_root+'/finley/lib') 22 sys.path.append(esys_root+'/escript/lib') 23 sys.path.append(esys_root+'/escript/py_src') 24 25 from escript import * 26 from util import * 27 from linearPDEs import * 28 from numarray import * 29 import finley as pdelib 30 31 32 # these values are currently fixed: 33 len_x0=1. 34 alpha=0.1 35 36 #############################################################################################################3 37 def solveVector(numDim, totalNumElem, len_x0, alpha, solver_method): 38 39 print "Vector solver:" 40 recDim=array([len_x0,1.,1.]) 41 # Define Computational Domain 42 numElem=int((totalNumElem/(len_x0*1.))**(1./numDim)) 43 elemDim = array([int(len_x0*numElem), numElem, numElem],Int) 44 45 # Set Mesh 46 if (numDim == 2): 47 mesh = pdelib.Rectangle(elemDim[0], elemDim[1], 2, \ 48 l0 = len_x0, l1 = 1.) 49 totElem=elemDim[0]*elemDim[1] 50 elif (numDim == 3): 51 mesh = pdelib.Brick(elemDim[0], elemDim[1], elemDim[2], 2, \ 52 l0 = len_x0, l1 = 1., l2 = 1.) 53 totElem=elemDim[0]*elemDim[1]*elemDim[2] 54 55 print " length of domain: ",recDim[:numDim] 56 print " requested elements: ",totalNumElem 57 print " num elements: ",totElem 58 # Set Mesh Descriptors 59 meshDim = mesh.getDim() 60 contfunc = ContinuousFunction(mesh) 61 func = Function(mesh) 62 x = contfunc.getX() 63 64 # Set Boundary Mask / pdelib Template "q" Parameter Vector 65 bndryMask = Vector(value = 0, what = contfunc) 66 for i in range(meshDim): 67 bndryMask += (x[i].whereZero() + (x[i]-recDim[i]).whereZero()) \ 68 * ones((numDim,)) 69 70 # Set True Solution / pdelib Template "r" Parameter Vector 71 u = Vector(value = 0, what = contfunc) 72 for i in range(meshDim): 73 for j in range(meshDim - 1): 74 u[i] += x[(i + j + 1) % meshDim]**2 75 # Set pdelib Template "A" Parameter Tensor 76 A = Tensor4(value = 0, what = func) 77 for i in range(meshDim): 78 for j in range(meshDim): 79 A[i,i,j,j] += 1. 80 A[i,j,j,i] += alpha 81 A[i,j,i,j] += alpha 82 83 # Build the pdelib System Matrix and RHS 84 mypde=LinearPDE(domain=mesh, \ 85 A = A, Y = - 2 * alpha * (meshDim - 1)*ones(meshDim), q = bndryMask, r = u) 86 mypde.setSolverMethod(solver_method) 87 88 # Solve for Approximate Solution 89 u_approx = mypde.getSolution(iter_max=10000) 90 91 # Report Results 92 error=Lsup(u - u_approx)/Lsup(u) 93 print " error L^sup Norm : ", error 94 print " residual L^sup Norm : ", Lsup(mypde.getResidual(u_approx)) 95 96 return error 97 98 ################################################################################################################# 99 100 def solveScalar(numDim, totalNumElem, len_x0, alpha, solver_method): 101 102 print "Scalar solver:" 103 recDim=array([len_x0,1.,1.]) 104 # Define Computational Domain 105 numElem=int((totalNumElem/(len_x0*1.))**(1./numDim)) 106 elemDim = array([int(len_x0*numElem), numElem, numElem],Int) 107 # Set Mesh 108 if (numDim == 2): 109 mesh = pdelib.Rectangle(elemDim[0], elemDim[1], 2, \ 110 l0 = len_x0, l1 = 1.) 111 totElem=elemDim[0]*elemDim[1] 112 elif (numDim == 3): 113 mesh = pdelib.Brick(elemDim[0], elemDim[1], elemDim[2], 2, \ 114 l0 = len_x0, l1 = 1., l2 = 1.) 115 totElem=elemDim[0]*elemDim[1]*elemDim[2] 116 117 print " length of domain: ",recDim[:numDim] 118 print " requested elements: ",totalNumElem 119 print " num elements: ",totElem 120 121 # Set Mesh Descriptors 122 meshDim = mesh.getDim() 123 contfunc = ContinuousFunction(mesh) 124 func = Function(mesh) 125 x = contfunc.getX() 126 127 # Set Boundary Mask / pdelib Template "q" Parameter Vector 128 bndryMask = Scalar(value = 0, what = contfunc) 129 for i in range(meshDim): 130 bndryMask += (x[i].whereZero() + (x[i]-recDim[i]).whereZero()) * 1.0 131 132 # Set True Solution / pdelib Template "r" Parameter Vector 133 u = Scalar(value = 0, what = contfunc) 134 for j in range(meshDim): 135 u += x[j] * x[j] 136 137 # Build the pdelib System Matrix and RHS 138 mypde=LinearPDE(domain=mesh, \ 139 A = identity(numDim), D = alpha, Y = alpha * u - 2 * meshDim, q = bndryMask, r = u) 140 mypde.setSolverMethod(solver_method) 141 142 # Solve for Approximate Solution 143 u_approx = mypde.getSolution(iter_max=10000) 144 145 # Report Results 146 error=Lsup(u - u_approx)/Lsup(u) 147 print " error L^sup Norm : ", error 148 print " residual L^sup Norm : ", Lsup(mypde.getResidual(u_approx)) 149 150 return error 151 152 ####################################################################################### 153 154 155 print "Test is started:" 156 print "----------------" 157 error=0. 158 for numDim in [2,3]: 159 for totalNumElem in [100, 200, 400, 800, 1600, 3200, 6400, 12800, 25600, 51200, 102400]: 160 for problem in [solveScalar,solveVector]: 161 # for solver_method in [ LinearPDE.PRES20, LinearPDE.PCG, LinearPDE.DIRECT, LinearPDE.BICGSTAB]: 162 for solver_method in [ LinearPDE.PRES20, LinearPDE.PCG, LinearPDE.BICGSTAB]: 163 error=max([problem(numDim, totalNumElem, len_x0, alpha, solver_method),error]) 164 print "----------------" 165 print "maximum error over all tests is ",error 166 print "----------------"

## Properties

Name Value
svn:eol-style native
svn:keywords Author Date Id Revision