C*************************************************************** C Test program for subroutine AEBEMA by Stephen Kirkup C*************************************************************** C C Copyright 1998- Stephen Kirkup C John Tyndall Nuclear Research Institute C School of Computing Engineering and Physical Sciences C University of Central Lancashire - www.uclan.ac.uk C Westlakes Campus C Samuel Lindow Building C West Lakes Science and Technology Park C Whitehaven C Cumbria CA24 3JY C United Kingdom C smkirkup@uclan.ac.uk C C This open source code can be found at C www.boundary-element-method.com/fortran/AEBEMA_T.FOR C C Issued under the GNU General Public License 2007, see gpl.txt C C Part of the the author's open source BEM packages. C All codes and manuals can be downloaded from C www.boundary-element-method.com C C*************************************************************** C This program is a test for the subroutine AEBEMA. The program computes C the solution to an acoustic/Helmholtz problem exterior to a sphere C by the boundary element method. C C Background C ---------- C C The Helmholtz problem arises when harmonic solutions of the wave C equation C 2 C __ 2 1 d {\Psi}(p,t) C \/ {\Psi}(p,t) - ---- --------------- = 0 C 2 2 C c d t C C are sought, where {\Psi}(p,t) is the scalar time-dependent velocity C potential. In the cases where {\Psi} is periodic, it may be C approximated as a set of frequency components that may be analysed C independently. For each frequency a component of the form C C {\phi}(p) exp(i {\omega} t) C C (where {\omega} = 2 * {\pi} * frequency) the wave equation can be C reduced to the Helmholtz equation C C __ 2 2 C \/ {\phi} + k {\phi} = 0 C C where k (the wavenumber) = {\omega}/c (c=speed of sound in the C medium). {\phi} is known as the velocity potential. C C For the exterior problem, the domain lies exterior to a closed C boundary S. The boundary condition may be Dirichlet, Robin or C Neumann. It is assumed to have the following general form C C {\alpha}(q) {\phi}(q) + {\beta}(q) v(q) = f(q) C C where {\phi}(q) is the velocity potential at the point q on S, v(q) C is the derivative of {\phi} with respect to the outward normal to S C at q and {\alpha}, {\beta} and f are complex-valued functions defined C on S. C C Subroutine AEBEMA accepts the wavenumber, a description of the C boundary of the domain and the position of the exterior points C where the solution ({\phi}) is sought, the boundary condition and C returns the solution ({\phi} and v) on S and the value of {\phi} C at the exterior points. C C The test problems C ----------------- C C In this test the domain is a sphere of diameter 1 (metre). The acoustic C medium is air (at 20 celcius and 1 atmosphere, c=344.0 (metres per C second), density {\rho}=1.205 (kilograms per cubic metre) and the C solution to the problem with a Dirichlet boundary condition C ({\alpha}=1, {\beta}=0) and with a Neumann boundary condition C ({\alpha}=0, beta=1) are sought. For both problems the frequency is C 400Hz (hence specifying k). C C In the R-z plane, the boundary conditions are specified through C taking the solution to be determined by C C {\phi} = sin(k z) , C C which is clearly a solution of the Helmholtz equation. C C C The boundary is described by a set of NS=32 elements of equal size, C so that each side comprises eight elements. The boundary solution C points are the centres of the elements. C The *s show the exterior points at which the solution is sought; C the points (0.025,0.025), (0.075,0.025), (0.025,0.075), C (0.075,0.075), and (0.05,0.05). C---------------------------------------------------------------------- C The PARAMETER statement C ----------------------- C There are four components in the PARAMETER statement. C integer MAXNS : The limit on the number of boundary elements. C integer MAXNV : The limit on the number of vertices. C integer MAXNPE : The limit on the number of exterior points. C External modules related to the package C --------------------------------------- C subroutine AEBEMA: Subroutine for solving the exterior Helmholtz C equation. (file AEBEMA.FOR contains AEBEMA and subordinate routines) C subroutine H2LC: Returns the individual discrete Helmholtz integral C operators. (file H2LC.FOR contains H2LC and subordinate routines) C subroutine CGLS: Solves a general linear system of equations. C (file CGLS.FOR contains CGSL and subordinate routines) C subroutine FNHANK: This computes Hankel functions of the first kind C and of order zero and one. (e.g. file FNHANK.FOR) C file GEOM2D.FOR contains the set of relevant geometric subroutines C The program PROGRAM AEBEMAT IMPLICIT NONE C VARIABLE DECLARATION C -------------------- C PARAMETERs for storing the limits on the dimension of arrays C Limit on the number of elements INTEGER MAXNS PARAMETER (MAXNS=32) C Limit on the number of vertices (equal to the number of elements) INTEGER MAXNV PARAMETER (MAXNV=MAXNS+1) C Limit on the number of test problems INTEGER MAXTEST PARAMETER (MAXTEST=200) C Limit on the number of points exterior to the boundary, where C acoustic properties are sought INTEGER MAXNPE PARAMETER (MAXNPE=1) C Constants C Real scalars: 0, 1, 2, pi REAL*8 ZERO,ONE,TWO,FOUR,PI C Complex scalars: (0,0), (1,0), (0,1) COMPLEX*16 CZERO,CONE,CIMAG C The reference pressure, used to convert units to decibels. REAL*8 PREREF C Properties of the acoustic medium C The speed of sound [standard unit: metres per second] REAL*8 CVAL(MAXTEST) C The density [standard unit: kilograms per cubic metre] REAL*8 RHOVAL(MAXTEST) C Wavenumber parameter for AEBEMA REAL*8 K C Angular frequency COMPLEX*16 OMEGA C Geometrical description of the boundary(ies) C Number of elements and counter INTEGER NS,IS C Number of collocation points (on S) and counter INTEGER NSP,ISP C Number of vetices and counter INTEGER NV,IV C Index of nodal coordinate for defining boundaries (standard unit is C metres) REAL*8 VERTEX(MAXNV,2) C The two nodes that define each element on the boundaries INTEGER SELV(MAXNS,2) C The points exterior to the boundary(ies) where the acoustic C properties are sought and the directional vectors at those points. C [Only necessary if an exterior solution is sought.] C Number of exterior points and counter INTEGER NPE,IPE C Coordinates of the exterior points REAL*8 PEXT (MAXNPE,2) C Number of test problems and counter INTEGER NTEST,ITEST C Data structures that contain the parameters that define the test C problems C The acoustic frequency for each test. FRVAL(i) is assigned the C acoustic frequency of the i-th test problem. REAL*8 FRVAL(MAXTEST) C The wavenumber for each test. KVAL(i) is assigned the wavenumber C of the i-th test problem. REAL*8 KVAL(MAXTEST) C The nature of the boundary condition is specified by assigning C values to the data structures SALVAL and SBEVAL. C SALVAL(i,j) is assigned the value of {\alpha} at the center of the C j-th element for the i-th test problem. COMPLEX*16 SALVAL(MAXTEST,MAXNS) C SBEVAL(i,j) is assigned the value of {\beta} at the center of the C j-th element for the i-th test problem. COMPLEX*16 SBEVAL(MAXTEST,MAXNS) C The actual boundary condition is specified by assigning values to C the data structure SFVAL. C SFVAL(i,j) is assigned the value of f at the center of the j-th C element for the i-th test problem. COMPLEX*16 SFVAL(MAXTEST,MAXNS) C The incident potential at the centres of the elements in each test C case COMPLEX*16 SPHIIN(MAXTEST,MAXNS) C The derivative of the incident potential at the centres of the C elements in each test COMPLEX*16 SVELIN(MAXTEST,MAXNS) C The incident potential at the exterior points COMPLEX*16 PPHIIN(MAXTEST,MAXNPE) C Data structures that are used to define each test problem in turn C and are input parameters to AEBEMA. C SALPHA(j) is assigned the value of {\alpha} at the centre of the C j-th element. COMPLEX*16 SALPHA(MAXNS) C SBETA(j) is assigned the value of {\beta} at the centre of the C j-th element. COMPLEX*16 SBETA(MAXNS) C SF(j) is assigned the value of f at the centre of the j-th element. COMPLEX*16 SF(MAXNS) C The incident potential at the centres of the elements in each test C case COMPLEX*16 SFFPHI(MAXNS) C The derivative of the incident potential at the centres of the C elements in each test COMPLEX*16 SFFVEL(MAXNS) C The incident potential at the interior points COMPLEX*16 PFFPHI(MAXNPE) C Validation and control parameters for AEBEMA C Switch for particular solution LOGICAL LSOL C Validation switch LOGICAL LVALID C The maximum absolute error in the parameters that describe the C geometry of the boundary. REAL*8 EGEOM C The parameter COMPLEX*16 MU C Output from subroutine AEBEMA C The velocity potential (phi - the solution) at the centres of the C elements COMPLEX*16 SPHI(MAXNS) C The normal derivative of the velocity potential at the centres of the C elements COMPLEX*16 SVEL(MAXNS) C The velocity potential (phi - the solution) at exterior points COMPLEX*16 PEPHI(MAXNPE) C Element areas REAL*8 ELAREA(MAXNS),AREACN C Computing radiation ratio REAL*8 POWER,BPOWER,RADRAT C Workspace for AEBEMA COMPLEX*16 WKSPC1(MAXNS,MAXNS) COMPLEX*16 WKSPC2(MAXNS,MAXNS) COMPLEX*16 WKSPC3(MAXNPE,MAXNS) COMPLEX*16 WKSPC4(MAXNPE,MAXNS) COMPLEX*16 WKSPC5(MAXNS) COMPLEX*16 WKSPC6(MAXNS) LOGICAL WKSPC7(MAXNS) C Acoustic properties. These data structures are appended after each C execution of AEBEMA and contain the numerical solution to the test C problems. C At the centres of the elements C Sound pressure [standard unit: newtons per sphere metre (or C pascals) and phase] COMPLEX*16 SPRESS(MAXTEST,MAXNS) C Velocity potential phi COMPLEX*16 SPHIVAL(MAXTEST,MAXNS) C Velocity (v) [standard unit: metres per second (and phase)] C At the centres of the elements COMPLEX*16 SV(MAXTEST,MAXNS) C Sound intensity [standard unit: watts per sphere metre] REAL*8 SINTY(MAXTEST,MAXNS) C At the exterior points C Sound pressure [standard unit: newtons per sphere metre (or C pascals) and phase] COMPLEX*16 IPRESS(MAXTEST,MAXNPE) C Velocity potential phi COMPLEX*16 IPHIVAL(MAXTEST,MAXNPE) C Counter through the x,y coordinates INTEGER ICOORD C Local storage of pressure, pressure/velocity COMPLEX*16 PRESSURE C The coordinates of the centres of the elements REAL*8 SELCNT(MAXNS,2) REAL*8 EPS C INITIALISATION C -------------- C Set constants ZERO=0.0D0 ONE=1.0D0 TWO=2.0D0 FOUR=4.0D0 PI=4.0D0*ATAN(ONE) CZERO=CMPLX(ZERO,ZERO) CONE=CMPLX(ONE,ZERO) CIMAG=CMPLX(ZERO,ONE) EPS=1.0E-10 C Reference for decibel scales PREREF=2.0D-05 C Describe the nodes and elements that make up the boundary C :The circle that generates the sphere is divided into NS=18 uniform C : elements. VERTEX and SELV are defined anti-clockwise around the C : boundary so that the normal to the boundary is assumed to be C : outward C :Set up nodes C : Set NS, the number of elements NS=18 C : Set NV, the number of vertices (equal to the number of elements) NV=NS+1 C : Set coordinates of the nodes DATA ((VERTEX(IV,ICOORD),ICOORD=1,2),IV=1,19) * / 0.000D0, 1.000D0, * 0.174D0, 0.985D0, * 0.342D0, 0.940D0, * 0.500D0, 0.866D0, * 0.643D0, 0.766D0, * 0.766D0, 0.643D0, * 0.866D0, 0.500D0, * 0.940D0, 0.342D0, * 0.985D0, 0.174D0, * 1.000D0, 0.000D0, * 0.985D0,-0.174D0, * 0.940D0,-0.342D0, * 0.866D0,-0.500D0, * 0.766D0,-0.643D0, * 0.643D0,-0.766D0, * 0.500D0,-0.866D0, * 0.342D0,-0.940D0, * 0.174D0,-0.985D0, * 0.000D0,-1.000D0 / C :Describe the elements that make up the two boundarys C : Set NS, the number of elements NS=18 C : Set nodal indices that describe the elements of the boundarys. C : The indices refer to the nodes in VERTEX. The order of the C : nodes in SELV dictates that the normal is outward from the C : boundary into the acoustic domain. DATA ((SELV(IS,ICOORD),ICOORD=1,2),IS=1,18) * / 1, 2, * 2, 3, * 3, 4, * 4, 5, * 5, 6, * 6, 7, * 7, 8, * 8, 9, * 9, 10, * 10, 11, * 11, 12, * 12, 13, * 13, 14, * 14, 15, * 15, 16, * 16, 17, * 17, 18, * 18, 19 / C Set the centres of the elements, the collocation points DO IS=1,NS SELCNT(IS,1)=(VERTEX(SELV(IS,1),1) * +VERTEX(SELV(IS,2),1))/TWO SELCNT(IS,2)=(VERTEX(SELV(IS,1),2) * +VERTEX(SELV(IS,2),2))/TWO ELAREA(IS)=AREACN(VERTEX(SELV(IS,1),1),VERTEX(SELV(IS,1),2), * VERTEX(SELV(IS,2),1),VERTEX(SELV(IS,2),2)) END DO C Set the points in the acoustic domain where the acoustic properties C are sought, PEXT . NPE=0 C The number of points on the boundary is equal to the number of C elements NSP=NS C Set up test problems C :Set the number of test problems NTEST=10 DO 400 ITEST=1,NTEST C Properties of the acoustic medium. C the propagation velocity C and RHO the density of the acoustic medium. C>0, RHO>0 C :Acoustic medium is air at 20 celcius and 1 atmosphere. C [C in metres per second, RHO in kilograms per cubic metre.] CVAL(ITEST)=344.0D0 RHOVAL(ITEST)=1.205D0 C :Set acoustic frequency value (hertz) in FRVAL FRVAL(ITEST)=ITEST*100.0D0 C : Set the wavenumber in KVAL KVAL(ITEST)=TWO*PI*FRVAL(ITEST)/CVAL(ITEST) C :Set nature of the boundary condition by prescribing the values of C the boundary functions SALVAL and SBEVAL at the collocation points C :In this case a Dirichlet (phi-valued) boundary condition K=KVAL(ITEST) DO 160 ISP=1,NSP SALVAL(ITEST,ISP)=CZERO SBEVAL(ITEST,ISP)=CONE SFVAL(ITEST,ISP)=1.0D0 160 CONTINUE DO 180 ISP=1,NSP SPHIIN(ITEST,ISP)=CZERO SVELIN(ITEST,ISP)=CZERO 180 CONTINUE DO 190 IPE=1,NPE PPHIIN(ITEST,IPE)=CZERO 190 CONTINUE C Set up validation and control parameters C :Switch for particular solution LSOL=.TRUE. C :Switch on the validation of AEBEMA LVALID=.TRUE. C :Set EGEOM EGEOM=1.0D-6 400 CONTINUE C Loop(ITEST) through the test problems DO 500 ITEST=1,NTEST C Set OMEGA, the angular frequency omega and K, the wavenumber K=KVAL(ITEST) OMEGA=2.0D0*PI*FRVAL(ITEST) C Set up particular alpha and beta functions for this wavenumber C and type of boundary condition DO 640 ISP=1,NSP SALPHA(ISP)=SALVAL(ITEST,ISP) SBETA(ISP)=SBEVAL(ITEST,ISP) SF(ISP)=SFVAL(ITEST,ISP) SFFPHI(ISP)=SPHIIN(ITEST,ISP) SFFVEL(ISP)=SVELIN(ITEST,ISP) 640 CONTINUE DO 650 IPE=1,NPE PFFPHI(IPE)=PPHIIN(ITEST,IPE) 650 CONTINUE C Set the parameter MU=CIMAG/(K+ONE) CALL AEBEMA(K, * MAXNV,NV,VERTEX,MAXNS,NS,SELV, * MAXNPE,NPE,PEXT , * SALPHA,SBETA,SF,SFFPHI,SFFVEL,PFFPHI, * LSOL,LVALID,EGEOM,MU, * SPHI,SVEL,PEPHI, * WKSPC1,WKSPC2,WKSPC3,WKSPC4, * WKSPC5,WKSPC6,WKSPC7) C Compute the sound pressure at the exterior points. Also compute C the velocity and intensity at the points for each type of boundary C condition and each related input function f and at each point. DO 690 ISP=1,NSP SPHIVAL(ITEST,ISP)=SPHI(ISP) PRESSURE=CIMAG*RHOVAL(ITEST)*OMEGA*SPHI(ISP) SPRESS(ITEST,ISP)=PRESSURE SV(ITEST,ISP)=SVEL(ISP) SINTY(ITEST,ISP)= * DBLE(CONJG(PRESSURE)*SVEL(ISP))/TWO 690 CONTINUE DO 695 ISP=1,NPE IPHIVAL(ITEST,ISP)=PEPHI(ISP) IPRESS(ITEST,ISP)=CIMAG*RHOVAL(ITEST)*OMEGA*PEPHI(ISP) 695 CONTINUE C Close loop(ITEST) through the test problems 500 CONTINUE C Output the solutions C Open file for the output data OPEN(UNIT=20,FILE='AEBEMA.OUT') C Formats for output 999 FORMAT(3F16.4) WRITE(20,*) 'AEBEMA: Computed solution to the acoustic properties' WRITE(20,*) WRITE(20,*) 'The radiation ratio of a pulsating sphere' WRITE(20,*) WRITE(20,*) WRITE(20,*) ' Frequency Wavenumber '// *'Radiation ratio' C Loop(ITEST) through the test problems. DO 2000 ITEST=1,NTEST C Output the acoustic frequency C Loop(ISP) through the points on the boundary POWER=0.0D0 BPOWER=0.0D0 DO 2030 ISP=1,NSP C Output the sound pressure, velocity and intensity at each point POWER=POWER+SINTY(ITEST,ISP)*ELAREA(ISP) BPOWER=BPOWER+RHOVAL(ITEST)*CVAL(ITEST)* * ABS(SV(ITEST,ISP))**2*ELAREA(ISP)/2 2030 CONTINUE RADRAT=POWER/BPOWER WRITE(20,999) FRVAL(ITEST),KVAL(ITEST),RADRAT 2000 CONTINUE CLOSE(20) END REAL*8 FUNCTION AREACN(R1,Z1,R2,Z2) REAL*8 R1,Z1,R2,Z2 REAL*8 A,B REAL*8 PI PI=4.0D0*ATAN(1.0D0) IF (ABS(Z2-Z1).GT.0.000001) THEN A=R1-(R2-R1)*Z1/(Z2-Z1) B=(R2-R1)/(Z2-Z1) AREACN= * ABS(2.0D0*PI*SQRT(1.0D0+B*B)*(Z2-Z1)*(A+B*(Z1+Z2)/2.0D0)) ELSE AREA=PI*ABS(R1*R1-R2*R2) END IF END