Contents of /branches/diaplayground/helmholtz.py

Revision 4944 - (show annotations)
Thu May 15 07:00:46 2014 UTC (5 years, 4 months ago) by caltinay
File MIME type: text/x-python
File size: 1255 byte(s)
```blocksize >1 lookin' good...

```
 1 from esys.escript import * 2 from esys.escript.linearPDEs import LinearPDE,SolverOptions 3 from esys.ripley import Rectangle 4 from time import time 5 6 BLOCKSIZE=4 7 8 dom = Rectangle(l0=1.,l1=1.,n0=2, n1=2) 9 x = dom.getX() 10 n = dom.getNormal() 11 12 pde = LinearPDE(dom, numEquations=BLOCKSIZE, numSolutions=BLOCKSIZE) 13 A = pde.createCoefficient("A") 14 D = pde.createCoefficient("D") 15 Y = pde.createCoefficient("Y") 16 d = pde.createCoefficient("d") 17 y = pde.createCoefficient("y") 18 if BLOCKSIZE == 1: 19 omega=0.1 20 A[:,:]=kronecker(dom) 21 D=omega 22 Y=omega*x[0] 23 d=10 24 y=n[0]+10*x[0] 25 else: 26 for i in range(BLOCKSIZE): 27 omega=0.1*i 28 A[i,:,i,:]=kronecker(dom) 29 D[i,:]=[omega]*BLOCKSIZE 30 Y[i]=omega*x[0] 31 d[i,:]=[10]*BLOCKSIZE 32 y[i]=n[0]+10*x[0] 33 34 pde.setValue(A=A, D=D, Y=Y, d=d, y=y) 35 pde.getSolverOptions().setSolverMethod(SolverOptions.PCG) 36 pde.getSolverOptions().setPreconditioner(SolverOptions.NO_PRECONDITIONER) 37 pde.getSolverOptions().setVerbosityOn() 38 #pde.setDebugOn() 39 rhs=pde.getRightHandSide() 40 #saveDataCSV('helmholtz_rhs.csv',rhs=rhs) 41 pde.getSystem()[0].saveMM('/tmp/helmholtzripley.mtx') 42 t0=time() 43 x = pde.getSolution() 44 t1=time() 45 print "Paso Solver Time: ", t1-t0 46 print "Solution: %s..%s"%(inf(x),sup(x)) 47 print x 48