# Annotation of /trunk/esys2/finley/test/python/SolveTest.py

Revision 126 - (hide annotations)
Fri Jul 22 03:53:08 2005 UTC (15 years, 8 months ago) by jgs
File MIME type: text/x-python
File size: 6274 byte(s)
```Merge of development branch back to main trunk on 2005-07-22

```
 1 jgs 102 # \$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 jgs 104 from esys.escript import * 18 from esys.linearPDEs import * 19 import esys.finley as pdelib 20 jgs 117 from time import time 21 jgs 104 22 jgs 102 from numarray import * 23 24 # these values are currently fixed: 25 len_x0=1. 26 jgs 117 alpha=10. 27 tm=0 28 jgs 102 29 #############################################################################################################3 30 jgs 115 def solveVector(numDim, totalNumElem, len_x0, alpha, solver_method,prec): 31 jgs 117 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 jgs 102 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 += (x[i].whereZero() + (x[i]-recDim[i]).whereZero()) \ 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 jgs 104 mypde=LinearPDE(mesh) 83 mypde.setValue(A = A, Y = - 2 * alpha * (meshDim - 1)*ones(meshDim), q = bndryMask, r = u) 84 jgs 117 mypde.setSolverMethod(solver_method_id) 85 jgs 123 # mypde.getOperator().saveMM("g.mm") 86 jgs 102 87 # Solve for Approximate Solution 88 jgs 117 tm=time() 89 u_approx = mypde.getSolution(preconditioner=prec_id,iter_max=10000) 90 tm=time()-tm 91 jgs 102 92 # Report Results 93 error=Lsup(u - u_approx)/Lsup(u) 94 jgs 117 print "@@ Vector %d : %d : %s(%s): error L^sup Norm : %e, time %e"%(mypde.getDim(),totElem,solver_method,prec,error,tm) 95 jgs 102 96 return error 97 98 ################################################################################################################# 99 100 jgs 115 def solveScalar(numDim, totalNumElem, len_x0, alpha, solver_method,prec): 101 jgs 102 102 jgs 117 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 jgs 102 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 += (x[i].whereZero() + (x[i]-recDim[i]).whereZero()) * 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 jgs 104 mypde=LinearPDE(mesh) 144 mypde.setValue(A = identity(numDim), D = alpha, Y = alpha * u - 2 * meshDim, q = bndryMask, r = u) 145 jgs 117 mypde.setSolverMethod(solver_method_id) 146 jgs 102 147 # Solve for Approximate Solution 148 jgs 117 tm=time() 149 u_approx = mypde.getSolution(preconditioner=prec_id,iter_max=10000) 150 tm=time()-tm 151 jgs 102 152 # Report Results 153 error=Lsup(u - u_approx)/Lsup(u) 154 jgs 117 print "@@ Scalar %d : %d : %s(%s): error L^sup Norm : %e, time %e"%(mypde.getDim(),totElem,solver_method,prec,error,tm) 155 jgs 102 156 return error 157 158 ####################################################################################### 159 160 161 print "Test is started:" 162 print "----------------" 163 error=0. 164 jgs 121 for numDim in [2, 3]: 165 jgs 123 # for totalNumElem in [51200]: 166 jgs 117 for totalNumElem in [100, 200, 400, 800, 1600, 3200, 6400, 12800, 25600, 51200, 102400,204800]: 167 for problem in [solveScalar,solveVector]: 168 jgs 123 #for problem in [solveVector]: 169 jgs 126 # error=max([problem(numDim, totalNumElem, len_x0, alpha,"PCG",""),error]) 170 jgs 122 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 jgs 102 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