# Contents of /trunk/dudley/test/python/RecTest.py

Revision 5706 - (show annotations)
Mon Jun 29 03:41:36 2015 UTC (4 years, 2 months ago) by sshaw
File MIME type: text/x-python
File size: 3851 byte(s)
all python files now force use of python3 prints and division syntax to stop sneaky errors appearing in py3 environs
 1 2 ############################################################################## 3 # 4 # Copyright (c) 2003-2015 by The University of Queensland 5 6 # 7 # Primary Business: Queensland, Australia 8 # Licensed under the Open Software License version 3.0 9 10 # 11 # Development until 2012 by Earth Systems Science Computational Center (ESSCC) 12 # Development 2012-2013 by School of Earth Sciences 13 # Development from 2014 by Centre for Geoscience Computing (GeoComp) 14 # 15 ############################################################################## 16 17 from __future__ import print_function, division 18 19 __copyright__="""Copyright (c) 2003-2015 by The University of Queensland 20 http://www.uq.edu.au 21 Primary Business: Queensland, Australia""" 22 __license__="""Licensed under the Open Software License version 3.0 23 24 __url__= 25 26 import sys 27 import os 28 import esys.escriptcore.utestselect as unittest 29 30 from esys.escript import * 31 from esys.escript.linearPDEs import * 32 from esys import dudley 33 34 import numpy 35 36 Pi=numpy.pi 37 numElements=15 38 sml=0.1 39 40 # 41 # this test the assemblage of problems with periodic boundary conditions: 42 # 43 44 # 45 # test solution is u_ex=sin(2*Pi*n*x0)*...*sin(2*Pi*n*x_dim) 46 # 47 def TheTest(msh,constraints,reduce): 48 n=ContinuousFunction(msh) 49 x=n.getX() 50 # set a test problem: 51 u_ex=Scalar(1,what=n) 52 for i in range(msh.getDim()): 53 u_ex*=sin(2*Pi*x[i]) 54 mypde=LinearPDE(msh) 55 mypde.setSymmetryOn() 56 mypde.setDebugOn() 57 mypde.setReducedOrderTo(reduce) 58 mypde.setValue(A=numpy.identity(msh.getDim()),D=sml,Y=(sml+4*Pi**2*msh.getDim())*u_ex,q=constraints,r=u_ex) 59 return Lsup(mypde.getSolution()-u_ex)/Lsup(u_ex) 60 61 max_error=0 62 max_text="none" 63 for onElements in [False,True]: 64 if onElements==True: 65 onElmtext=", with elements on faces" 66 else: 67 onElmtext="" 68 for order in [1,2]: 69 for reduce in [False,True]: 70 if reduce==True: 71 redtext=",reduced" 72 else: 73 redtext="" 74 for dim in [2,3]: 75 if dim==2: 76 for i0 in [True,True]: 77 for i1 in [True,True]: 78 msh=dudley.Rectangle(numElements,numElements,order,periodic0=i0,periodic1=i1,useElementsOnFace=onElements) 79 n=ContinuousFunction(msh) 80 x=n.getX() 81 c=Scalar(0,what=n) 82 if i0==False: 83 c+=whereZero(x[0])+whereZero(x[0]-1.) 84 if i1==False: 85 c+=whereZero(x[1])+whereZero(x[1]-1.) 86 error=TheTest(msh,c,reduce) 87 text="Rectangle order = %d%s%s, periodic0= %d, periodic1= %d: error= %f"%(order,redtext,onElmtext,i0,i1,error) 88 print("@@ ",text) 89 if error>max_error: 90 max_error=error 91 max_text=text 92 elif dim==3: 93 for i0 in [True,False]: 94 for i1 in [True,False]: 95 for i2 in [True,False]: 96 msh=dudley.Brick(numElements,numElements,numElements,order,periodic0=i0,periodic1=i1,periodic2=i2,useElementsOnFace=onElements) 97 n=ContinuousFunction(msh) 98 x=n.getX() 99 c=Scalar(0,what=n) 100 if i0==False: 101 c+=whereZero(x[0])+whereZero(x[0]-1.) 102 if i1==False: 103 c+=whereZero(x[1])+whereZero(x[1]-1.) 104 if i2==False: 105 c+=whereZero(x[2])+whereZero(x[2]-1.) 106 error=TheTest(msh,c,reduce) 107 text="Brick order = %d%s%s, periodic0= %d, periodic1= %d, periodic2= %d: error= %f"%(order,redtext,onElmtext,i0,i1,i2,error) 108 print("@@ ",text) 109 if error>max_error: 110 max_error=error 111 max_text=text 112 113 print("@@@@ maximum error for :",max_text)

## Properties

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

 ViewVC Help Powered by ViewVC 1.1.26