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

Revision 147 - (show annotations)
Fri Aug 12 01:45:47 2005 UTC (14 years, 7 months ago) by jgs
Original Path: trunk/esys2/finley/test/python/SolveTest.py
File MIME type: text/x-python
File size: 6372 byte(s)
```erge of development branch dev-02 back to main trunk on 2005-08-12

```
 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 # from esys.escript import * 18 # from esys.linearPDEs import * 19 # import esys.finley as pdelib 20 from escript.escript import * 21 from escript.linearPDEs import * 22 import finley.finley as pdelib 23 from time import time 24 25 from numarray import * 26 27 # these values are currently fixed: 28 len_x0=1. 29 alpha=10. 30 tm=0 31 32 #############################################################################################################3 33 def solveVector(numDim, totalNumElem, len_x0, alpha, solver_method,prec): 34 35 if prec=="": 36 prec_id=0 37 else: 38 prec_id=eval("LinearPDE.%s"%prec) 39 solver_method_id=eval("LinearPDE.%s"%solver_method) 40 print "Vector solver:" 41 recDim=array([len_x0,1.,1.]) 42 # Define Computational Domain 43 numElem=int((totalNumElem/(len_x0*1.))**(1./numDim)) 44 elemDim = array([int(len_x0*numElem), numElem, numElem],Int) 45 46 # Set Mesh 47 if (numDim == 2): 48 mesh = pdelib.Rectangle(elemDim[0], elemDim[1], 2, \ 49 l0 = len_x0, l1 = 1.) 50 totElem=elemDim[0]*elemDim[1] 51 elif (numDim == 3): 52 mesh = pdelib.Brick(elemDim[0], elemDim[1], elemDim[2], 2, \ 53 l0 = len_x0, l1 = 1., l2 = 1.) 54 totElem=elemDim[0]*elemDim[1]*elemDim[2] 55 56 print " length of domain: ",recDim[:numDim] 57 print " requested elements: ",totalNumElem 58 print " num elements: ",totElem 59 # Set Mesh Descriptors 60 meshDim = mesh.getDim() 61 contfunc = ContinuousFunction(mesh) 62 func = Function(mesh) 63 x = contfunc.getX() 64 65 # Set Boundary Mask / pdelib Template "q" Parameter Vector 66 bndryMask = Vector(value = 0, what = contfunc) 67 for i in range(meshDim): 68 bndryMask += (x[i].whereZero() + (x[i]-recDim[i]).whereZero()) \ 69 * ones((numDim,)) 70 71 # Set True Solution / pdelib Template "r" Parameter Vector 72 u = Vector(value = 0, what = contfunc) 73 for i in range(meshDim): 74 for j in range(meshDim - 1): 75 u[i] += x[(i + j + 1) % meshDim]**2 76 # Set pdelib Template "A" Parameter Tensor 77 A = Tensor4(value = 0, what = func) 78 for i in range(meshDim): 79 for j in range(meshDim): 80 A[i,i,j,j] += 1. 81 A[i,j,j,i] += alpha 82 A[i,j,i,j] += alpha 83 84 # Build the pdelib System Matrix and RHS 85 mypde=LinearPDE(mesh) 86 mypde.setValue(A = A, Y = - 2 * alpha * (meshDim - 1)*ones(meshDim), q = bndryMask, r = u) 87 mypde.setSolverMethod(solver_method_id) 88 # mypde.getOperator().saveMM("g.mm") 89 90 # Solve for Approximate Solution 91 tm=time() 92 u_approx = mypde.getSolution(preconditioner=prec_id,iter_max=10000) 93 tm=time()-tm 94 95 # Report Results 96 error=Lsup(u - u_approx)/Lsup(u) 97 print "@@ Vector %d : %d : %s(%s): error L^sup Norm : %e, time %e"%(mypde.getDim(),totElem,solver_method,prec,error,tm) 98 99 return error 100 101 ################################################################################################################# 102 103 def solveScalar(numDim, totalNumElem, len_x0, alpha, solver_method,prec): 104 105 if prec=="": 106 prec_id=0 107 else: 108 prec_id=eval("LinearPDE.%s"%prec) 109 solver_method_id=eval("LinearPDE.%s"%solver_method) 110 print "Scalar solver:" 111 recDim=array([len_x0,1.,1.]) 112 # Define Computational Domain 113 numElem=int((totalNumElem/(len_x0*1.))**(1./numDim)) 114 elemDim = array([int(len_x0*numElem), numElem, numElem],Int) 115 # Set Mesh 116 if (numDim == 2): 117 mesh = pdelib.Rectangle(elemDim[0], elemDim[1], 2, \ 118 l0 = len_x0, l1 = 1.) 119 totElem=elemDim[0]*elemDim[1] 120 elif (numDim == 3): 121 mesh = pdelib.Brick(elemDim[0], elemDim[1], elemDim[2], 2, \ 122 l0 = len_x0, l1 = 1., l2 = 1.) 123 totElem=elemDim[0]*elemDim[1]*elemDim[2] 124 125 print " length of domain: ",recDim[:numDim] 126 print " requested elements: ",totalNumElem 127 print " num elements: ",totElem 128 129 # Set Mesh Descriptors 130 meshDim = mesh.getDim() 131 contfunc = ContinuousFunction(mesh) 132 func = Function(mesh) 133 x = contfunc.getX() 134 135 # Set Boundary Mask / pdelib Template "q" Parameter Vector 136 bndryMask = Scalar(value = 0, what = contfunc) 137 for i in range(meshDim): 138 bndryMask += (x[i].whereZero() + (x[i]-recDim[i]).whereZero()) * 1.0 139 140 # Set True Solution / pdelib Template "r" Parameter Vector 141 u = Scalar(value = 0, what = contfunc) 142 for j in range(meshDim): 143 u += x[j] * x[j] 144 145 # Build the pdelib System Matrix and RHS 146 mypde=LinearPDE(mesh) 147 mypde.setValue(A = identity(numDim), D = alpha, Y = alpha * u - 2 * meshDim, q = bndryMask, r = u) 148 mypde.setSolverMethod(solver_method_id) 149 150 # Solve for Approximate Solution 151 tm=time() 152 u_approx = mypde.getSolution(preconditioner=prec_id,iter_max=10000) 153 tm=time()-tm 154 155 # Report Results 156 error=Lsup(u - u_approx)/Lsup(u) 157 print "@@ Scalar %d : %d : %s(%s): error L^sup Norm : %e, time %e"%(mypde.getDim(),totElem,solver_method,prec,error,tm) 158 159 return error 160 161 ####################################################################################### 162 163 164 print "Test is started:" 165 print "----------------" 166 error=0. 167 for numDim in [2, 3]: 168 # for totalNumElem in [51200]: 169 for totalNumElem in [100, 200, 400, 800, 1600, 3200, 6400, 12800, 25600, 51200, 102400,204800]: 170 for problem in [solveScalar,solveVector]: 171 #for problem in [solveVector]: 172 error=max([problem(numDim, totalNumElem, len_x0, alpha,"PCG",""),error]) 173 error=max([problem(numDim, totalNumElem, len_x0, alpha,"DIRECT",""),error]) 174 #if totalNumElem*2**numDim*numDim< 200000: error=max([problem(numDim, totalNumElem, len_x0, alpha,"DIRECT",""),error]) 175 # for solver_method in [ "PCG" ]: 176 # for prec in [ "JACOBI", "ILU0" ]: 177 # error=max([problem(numDim, totalNumElem, len_x0, alpha, solver_method,prec),error]) 178 print "----------------" 179 print "maximum error over all tests is ",error 180 print "----------------"

## Properties

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