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

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

## Properties

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