Home

Context Navigation

source: palm/trunk/SOURCE/tridia_solver.f90 @ 1220

Last change on this file since 1220 was 1217, checked in by raasch, 11 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 18.9 KB

Rev	Line
[1212]	1	MODULE tridia_solver
	2
	3	!--------------------------------------------------------------------------------!
	4	! This file is part of PALM.
	5	!
	6	! PALM is free software: you can redistribute it and/or modify it under the terms
	7	! of the GNU General Public License as published by the Free Software Foundation,
	8	! either version 3 of the License, or (at your option) any later version.
	9	!
	10	! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
	11	! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
	12	! A PARTICULAR PURPOSE. See the GNU General Public License for more details.
	13	!
	14	! You should have received a copy of the GNU General Public License along with
	15	! PALM. If not, see <http://www.gnu.org/licenses/>.
	16	!
	17	! Copyright 1997-2012 Leibniz University Hannover
	18	!--------------------------------------------------------------------------------!
	19	!
	20	! Current revisions:
	21	! ------------------
[1213]	22	!
[1217]	23	!
[1213]	24	! Former revisions:
	25	! -----------------
	26	! $Id: tridia_solver.f90 1217 2013-08-26 11:30:05Z heinze $
	27	!
[1217]	28	! 1216 2013-08-26 09:31:42Z raasch
	29	! +tridia_substi_overlap for handling overlapping fft / transposition
	30	!
[1213]	31	! 1212 2013-08-15 08:46:27Z raasch
[1212]	32	! Initial revision.
	33	! Routines have been moved to seperate module from former file poisfft to here.
	34	! The tridiagonal matrix coefficients of array tri are calculated only once at
	35	! the beginning, i.e. routine split is called within tridia_init.
	36	!
	37	!
	38	! Description:
	39	! ------------
	40	! solves the linear system of equations:
	41	!
	42	! -(4 pi^2(i^2/(dx^2nnx^2)+j^2/(dy^2nny^2))+
	43	! 1/(dzu(k)dzw(k))+1/(dzu(k-1)dzw(k)))*p(i,j,k)+
	44	! 1/(dzu(k)dzw(k))p(i,j,k+1)+1/(dzu(k-1)dzw(k))p(i,j,k-1)=d(i,j,k)
	45	!
	46	! by using the Thomas algorithm
	47	!------------------------------------------------------------------------------!
	48
	49	USE indices
	50	USE transpose_indices
	51
	52	IMPLICIT NONE
	53
	54	REAL, DIMENSION(:,:), ALLOCATABLE :: ddzuw
	55
	56	PRIVATE
	57
	58	INTERFACE tridia_substi
	59	MODULE PROCEDURE tridia_substi
	60	END INTERFACE tridia_substi
	61
[1216]	62	INTERFACE tridia_substi_overlap
	63	MODULE PROCEDURE tridia_substi_overlap
	64	END INTERFACE tridia_substi_overlap
[1212]	65
[1216]	66	PUBLIC tridia_substi, tridia_substi_overlap, tridia_init, tridia_1dd
	67
[1212]	68	CONTAINS
	69
	70
	71	SUBROUTINE tridia_init
	72
	73	USE arrays_3d, ONLY: ddzu_pres, ddzw
	74
	75	IMPLICIT NONE
	76
	77	INTEGER :: k
	78
	79	ALLOCATE( ddzuw(0:nz-1,3) )
	80
	81	DO k = 0, nz-1
	82	ddzuw(k,1) = ddzu_pres(k+1) * ddzw(k+1)
	83	ddzuw(k,2) = ddzu_pres(k+2) * ddzw(k+1)
	84	ddzuw(k,3) = -1.0 * &
	85	( ddzu_pres(k+2) * ddzw(k+1) + ddzu_pres(k+1) * ddzw(k+1) )
	86	ENDDO
	87	!
	88	!-- Calculate constant coefficients of the tridiagonal matrix
	89	#if ! defined ( __check )
	90	CALL maketri
	91	CALL split
	92	#endif
	93
	94	END SUBROUTINE tridia_init
	95
	96
	97	SUBROUTINE maketri
	98
	99	!------------------------------------------------------------------------------!
	100	! Computes the i- and j-dependent component of the matrix
	101	!------------------------------------------------------------------------------!
	102
	103	USE arrays_3d, ONLY: tric
	104	USE constants
	105	USE control_parameters
	106	USE grid_variables
	107
	108	IMPLICIT NONE
	109
	110	INTEGER :: i, j, k, nnxh, nnyh
	111
	112	!$acc declare create( ll )
	113	REAL :: ll(nxl_z:nxr_z,nys_z:nyn_z)
	114
	115
	116	nnxh = ( nx + 1 ) / 2
	117	nnyh = ( ny + 1 ) / 2
	118
	119	!
	120	!-- Provide the constant coefficients of the tridiagonal matrix for solution
	121	!-- of the Poisson equation in Fourier space.
	122	!-- The coefficients are computed following the method of
	123	!-- Schmidt et al. (DFVLR-Mitteilung 84-15), which departs from Stephan
	124	!-- Siano's original version by discretizing the Poisson equation,
	125	!-- before it is Fourier-transformed.
	126
	127	!$acc kernels present( tric )
	128	!$acc loop vector( 32 )
	129	DO j = nys_z, nyn_z
	130	DO i = nxl_z, nxr_z
	131	IF ( j >= 0 .AND. j <= nnyh ) THEN
	132	IF ( i >= 0 .AND. i <= nnxh ) THEN
	133	ll(i,j) = 2.0 * ( 1.0 - COS( ( 2.0 * pi * i ) / &
	134	REAL( nx+1 ) ) ) / ( dx * dx ) + &
	135	2.0 * ( 1.0 - COS( ( 2.0 * pi * j ) / &
	136	REAL( ny+1 ) ) ) / ( dy * dy )
	137	ELSE
	138	ll(i,j) = 2.0 * ( 1.0 - COS( ( 2.0 * pi * ( nx+1-i ) ) / &
	139	REAL( nx+1 ) ) ) / ( dx * dx ) + &
	140	2.0 * ( 1.0 - COS( ( 2.0 * pi * j ) / &
	141	REAL( ny+1 ) ) ) / ( dy * dy )
	142	ENDIF
	143	ELSE
	144	IF ( i >= 0 .AND. i <= nnxh ) THEN
	145	ll(i,j) = 2.0 * ( 1.0 - COS( ( 2.0 * pi * i ) / &
	146	REAL( nx+1 ) ) ) / ( dx * dx ) + &
	147	2.0 * ( 1.0 - COS( ( 2.0 * pi * ( ny+1-j ) ) / &
	148	REAL( ny+1 ) ) ) / ( dy * dy )
	149	ELSE
	150	ll(i,j) = 2.0 * ( 1.0 - COS( ( 2.0 * pi * ( nx+1-i ) ) / &
	151	REAL( nx+1 ) ) ) / ( dx * dx ) + &
	152	2.0 * ( 1.0 - COS( ( 2.0 * pi * ( ny+1-j ) ) / &
	153	REAL( ny+1 ) ) ) / ( dy * dy )
	154	ENDIF
	155	ENDIF
	156	ENDDO
	157	ENDDO
	158
	159	!$acc loop
	160	DO k = 0, nz-1
	161	DO j = nys_z, nyn_z
	162	!$acc loop vector( 32 )
	163	DO i = nxl_z, nxr_z
	164	tric(i,j,k) = ddzuw(k,3) - ll(i,j)
	165	ENDDO
	166	ENDDO
	167	ENDDO
	168	!$acc end kernels
	169
	170	IF ( ibc_p_b == 1 ) THEN
	171	!$acc kernels present( tric )
	172	!$acc loop
	173	DO j = nys_z, nyn_z
	174	DO i = nxl_z, nxr_z
	175	tric(i,j,0) = tric(i,j,0) + ddzuw(0,1)
	176	ENDDO
	177	ENDDO
	178	!$acc end kernels
	179	ENDIF
	180	IF ( ibc_p_t == 1 ) THEN
	181	!$acc kernels present( tric )
	182	!$acc loop
	183	DO j = nys_z, nyn_z
	184	DO i = nxl_z, nxr_z
	185	tric(i,j,nz-1) = tric(i,j,nz-1) + ddzuw(nz-1,2)
	186	ENDDO
	187	ENDDO
	188	!$acc end kernels
	189	ENDIF
	190
	191	END SUBROUTINE maketri
	192
	193
	194	SUBROUTINE tridia_substi( ar )
	195
	196	!------------------------------------------------------------------------------!
	197	! Substitution (Forward and Backward) (Thomas algorithm)
	198	!------------------------------------------------------------------------------!
	199
	200	USE arrays_3d, ONLY: tri
	201	USE control_parameters
	202
	203	IMPLICIT NONE
	204
	205	INTEGER :: i, j, k
	206
	207	REAL :: ar(nxl_z:nxr_z,nys_z:nyn_z,1:nz)
	208
	209	!$acc declare create( ar1 )
	210	REAL, DIMENSION(nxl_z:nxr_z,nys_z:nyn_z,0:nz-1) :: ar1
	211
	212	!
	213	!-- Forward substitution
	214	DO k = 0, nz - 1
	215	!$acc kernels present( ar, tri )
	216	!$acc loop
	217	DO j = nys_z, nyn_z
	218	DO i = nxl_z, nxr_z
	219
	220	IF ( k == 0 ) THEN
	221	ar1(i,j,k) = ar(i,j,k+1)
	222	ELSE
	223	ar1(i,j,k) = ar(i,j,k+1) - tri(i,j,k,2) * ar1(i,j,k-1)
	224	ENDIF
	225
	226	ENDDO
	227	ENDDO
	228	!$acc end kernels
	229	ENDDO
	230
	231	!
	232	!-- Backward substitution
	233	!-- Note, the 1.0E-20 in the denominator is due to avoid divisions
	234	!-- by zero appearing if the pressure bc is set to neumann at the top of
	235	!-- the model domain.
	236	DO k = nz-1, 0, -1
	237	!$acc kernels present( ar, tri )
	238	!$acc loop
	239	DO j = nys_z, nyn_z
	240	DO i = nxl_z, nxr_z
	241
	242	IF ( k == nz-1 ) THEN
	243	ar(i,j,k+1) = ar1(i,j,k) / ( tri(i,j,k,1) + 1.0E-20 )
	244	ELSE
	245	ar(i,j,k+1) = ( ar1(i,j,k) - ddzuw(k,2) * ar(i,j,k+2) ) &
	246	/ tri(i,j,k,1)
	247	ENDIF
	248	ENDDO
	249	ENDDO
	250	!$acc end kernels
	251	ENDDO
	252
	253	!
	254	!-- Indices i=0, j=0 correspond to horizontally averaged pressure.
	255	!-- The respective values of ar should be zero at all k-levels if
	256	!-- acceleration of horizontally averaged vertical velocity is zero.
	257	IF ( ibc_p_b == 1 .AND. ibc_p_t == 1 ) THEN
	258	IF ( nys_z == 0 .AND. nxl_z == 0 ) THEN
	259	!$acc kernels loop present( ar )
	260	DO k = 1, nz
	261	ar(nxl_z,nys_z,k) = 0.0
	262	ENDDO
	263	ENDIF
	264	ENDIF
	265
	266	END SUBROUTINE tridia_substi
	267
	268
[1216]	269	SUBROUTINE tridia_substi_overlap( ar, jj )
	270
	271	!------------------------------------------------------------------------------!
	272	! Substitution (Forward and Backward) (Thomas algorithm)
	273	!------------------------------------------------------------------------------!
	274
	275	USE arrays_3d, ONLY: tri
	276	USE control_parameters
	277
	278	IMPLICIT NONE
	279
	280	INTEGER :: i, j, jj, k
	281
	282	REAL :: ar(nxl_z:nxr_z,nys_z:nyn_z,1:nz)
	283
	284	!$acc declare create( ar1 )
	285	REAL, DIMENSION(nxl_z:nxr_z,nys_z:nyn_z,0:nz-1) :: ar1
	286
	287	!
	288	!-- Forward substitution
	289	DO k = 0, nz - 1
	290	!$acc kernels present( ar, tri )
	291	!$acc loop
	292	DO j = nys_z, nyn_z
	293	DO i = nxl_z, nxr_z
	294
	295	IF ( k == 0 ) THEN
	296	ar1(i,j,k) = ar(i,j,k+1)
	297	ELSE
	298	ar1(i,j,k) = ar(i,j,k+1) - tri(i,jj,k,2) * ar1(i,j,k-1)
	299	ENDIF
	300
	301	ENDDO
	302	ENDDO
	303	!$acc end kernels
	304	ENDDO
	305
	306	!
	307	!-- Backward substitution
	308	!-- Note, the 1.0E-20 in the denominator is due to avoid divisions
	309	!-- by zero appearing if the pressure bc is set to neumann at the top of
	310	!-- the model domain.
	311	DO k = nz-1, 0, -1
	312	!$acc kernels present( ar, tri )
	313	!$acc loop
	314	DO j = nys_z, nyn_z
	315	DO i = nxl_z, nxr_z
	316
	317	IF ( k == nz-1 ) THEN
	318	ar(i,j,k+1) = ar1(i,j,k) / ( tri(i,jj,k,1) + 1.0E-20 )
	319	ELSE
	320	ar(i,j,k+1) = ( ar1(i,j,k) - ddzuw(k,2) * ar(i,j,k+2) ) &
	321	/ tri(i,jj,k,1)
	322	ENDIF
	323	ENDDO
	324	ENDDO
	325	!$acc end kernels
	326	ENDDO
	327
	328	!
	329	!-- Indices i=0, j=0 correspond to horizontally averaged pressure.
	330	!-- The respective values of ar should be zero at all k-levels if
	331	!-- acceleration of horizontally averaged vertical velocity is zero.
	332	IF ( ibc_p_b == 1 .AND. ibc_p_t == 1 ) THEN
	333	IF ( nys_z == 0 .AND. nxl_z == 0 ) THEN
	334	!$acc kernels loop present( ar )
	335	DO k = 1, nz
	336	ar(nxl_z,nys_z,k) = 0.0
	337	ENDDO
	338	ENDIF
	339	ENDIF
	340
	341	END SUBROUTINE tridia_substi_overlap
	342
	343
[1212]	344	SUBROUTINE split
	345
	346	!------------------------------------------------------------------------------!
	347	! Splitting of the tridiagonal matrix (Thomas algorithm)
	348	!------------------------------------------------------------------------------!
	349
	350	USE arrays_3d, ONLY: tri, tric
	351
	352	IMPLICIT NONE
	353
	354	INTEGER :: i, j, k
	355
	356	!
	357	!-- Splitting
	358	!$acc kernels present( tri, tric )
	359	!$acc loop
	360	DO j = nys_z, nyn_z
	361	!$acc loop vector( 32 )
	362	DO i = nxl_z, nxr_z
	363	tri(i,j,0,1) = tric(i,j,0)
	364	ENDDO
	365	ENDDO
	366	!$acc end kernels
	367
	368	DO k = 1, nz-1
	369	!$acc kernels present( tri, tric )
	370	!$acc loop
	371	DO j = nys_z, nyn_z
	372	!$acc loop vector( 32 )
	373	DO i = nxl_z, nxr_z
	374	tri(i,j,k,2) = ddzuw(k,1) / tri(i,j,k-1,1)
	375	tri(i,j,k,1) = tric(i,j,k) - ddzuw(k-1,2) * tri(i,j,k,2)
	376	ENDDO
	377	ENDDO
	378	!$acc end kernels
	379	ENDDO
	380
	381	END SUBROUTINE split
	382
	383
	384	SUBROUTINE tridia_1dd( ddx2, ddy2, nx, ny, j, ar, tri )
	385
	386	!------------------------------------------------------------------------------!
	387	! Solves the linear system of equations for a 1d-decomposition along x (see
	388	! tridia)
	389	!
	390	! Attention: when using the intel compilers older than 12.0, array tri must
	391	! be passed as an argument to the contained subroutines. Otherwise
	392	! addres faults will occur. This feature can be activated with
	393	! cpp-switch __intel11
	394	! On NEC, tri should not be passed (except for routine substi_1dd)
	395	! because this causes very bad performance.
	396	!------------------------------------------------------------------------------!
	397
	398	USE arrays_3d
	399	USE control_parameters
	400
	401	USE pegrid
	402
	403	IMPLICIT NONE
	404
	405	INTEGER :: i, j, k, nnyh, nx, ny, omp_get_thread_num, tn
	406
	407	REAL :: ddx2, ddy2
	408
	409	REAL, DIMENSION(0:nx,1:nz) :: ar
	410	REAL, DIMENSION(5,0:nx,0:nz-1) :: tri
	411
	412
	413	nnyh = ( ny + 1 ) / 2
	414
	415	!
	416	!-- Define constant elements of the tridiagonal matrix.
	417	!-- The compiler on SX6 does loop exchange. If 0:nx is a high power of 2,
	418	!-- the exchanged loops create bank conflicts. The following directive
	419	!-- prohibits loop exchange and the loops perform much better.
	420	! tn = omp_get_thread_num()
	421	! WRITE( 120+tn, * ) '+++ id=',myid,' nx=',nx,' thread=', omp_get_thread_num()
	422	! CALL local_flush( 120+tn )
	423	!CDIR NOLOOPCHG
	424	DO k = 0, nz-1
	425	DO i = 0,nx
	426	tri(2,i,k) = ddzu_pres(k+1) * ddzw(k+1)
	427	tri(3,i,k) = ddzu_pres(k+2) * ddzw(k+1)
	428	ENDDO
	429	ENDDO
	430	! WRITE( 120+tn, * ) '+++ id=',myid,' end of first tridia loop thread=', omp_get_thread_num()
	431	! CALL local_flush( 120+tn )
	432
	433	IF ( j <= nnyh ) THEN
	434	#if defined( __intel11 )
	435	CALL maketri_1dd( j, tri )
	436	#else
	437	CALL maketri_1dd( j )
	438	#endif
	439	ELSE
	440	#if defined( __intel11 )
	441	CALL maketri_1dd( ny+1-j, tri )
	442	#else
	443	CALL maketri_1dd( ny+1-j )
	444	#endif
	445	ENDIF
	446	#if defined( __intel11 )
	447	CALL split_1dd( tri )
	448	#else
	449	CALL split_1dd
	450	#endif
	451	CALL substi_1dd( ar, tri )
	452
	453	CONTAINS
	454
	455	#if defined( __intel11 )
	456	SUBROUTINE maketri_1dd( j, tri )
	457	#else
	458	SUBROUTINE maketri_1dd( j )
	459	#endif
	460
	461	!------------------------------------------------------------------------------!
	462	! computes the i- and j-dependent component of the matrix
	463	!------------------------------------------------------------------------------!
	464
	465	USE constants
	466
	467	IMPLICIT NONE
	468
	469	INTEGER :: i, j, k, nnxh
	470	REAL :: a, c
	471
	472	REAL, DIMENSION(0:nx) :: l
	473
	474	#if defined( __intel11 )
	475	REAL, DIMENSION(5,0:nx,0:nz-1) :: tri
	476	#endif
	477
	478
	479	nnxh = ( nx + 1 ) / 2
	480	!
	481	!-- Provide the tridiagonal matrix for solution of the Poisson equation in
	482	!-- Fourier space. The coefficients are computed following the method of
	483	!-- Schmidt et al. (DFVLR-Mitteilung 84-15), which departs from Stephan
	484	!-- Siano's original version by discretizing the Poisson equation,
	485	!-- before it is Fourier-transformed
	486	DO i = 0, nx
	487	IF ( i >= 0 .AND. i <= nnxh ) THEN
	488	l(i) = 2.0 * ( 1.0 - COS( ( 2.0 * pi * i ) / &
	489	REAL( nx+1 ) ) ) * ddx2 + &
	490	2.0 * ( 1.0 - COS( ( 2.0 * pi * j ) / &
	491	REAL( ny+1 ) ) ) * ddy2
	492	ELSE
	493	l(i) = 2.0 * ( 1.0 - COS( ( 2.0 * pi * ( nx+1-i ) ) / &
	494	REAL( nx+1 ) ) ) * ddx2 + &
	495	2.0 * ( 1.0 - COS( ( 2.0 * pi * j ) / &
	496	REAL( ny+1 ) ) ) * ddy2
	497	ENDIF
	498	ENDDO
	499
	500	DO k = 0, nz-1
	501	DO i = 0, nx
	502	a = -1.0 * ddzu_pres(k+2) * ddzw(k+1)
	503	c = -1.0 * ddzu_pres(k+1) * ddzw(k+1)
	504	tri(1,i,k) = a + c - l(i)
	505	ENDDO
	506	ENDDO
	507	IF ( ibc_p_b == 1 ) THEN
	508	DO i = 0, nx
	509	tri(1,i,0) = tri(1,i,0) + tri(2,i,0)
	510	ENDDO
	511	ENDIF
	512	IF ( ibc_p_t == 1 ) THEN
	513	DO i = 0, nx
	514	tri(1,i,nz-1) = tri(1,i,nz-1) + tri(3,i,nz-1)
	515	ENDDO
	516	ENDIF
	517
	518	END SUBROUTINE maketri_1dd
	519
	520
	521	#if defined( __intel11 )
	522	SUBROUTINE split_1dd( tri )
	523	#else
	524	SUBROUTINE split_1dd
	525	#endif
	526
	527	!------------------------------------------------------------------------------!
	528	! Splitting of the tridiagonal matrix (Thomas algorithm)
	529	!------------------------------------------------------------------------------!
	530
	531	IMPLICIT NONE
	532
	533	INTEGER :: i, k
	534
	535	#if defined( __intel11 )
	536	REAL, DIMENSION(5,0:nx,0:nz-1) :: tri
	537	#endif
	538
	539
	540	!
	541	!-- Splitting
	542	DO i = 0, nx
	543	tri(4,i,0) = tri(1,i,0)
	544	ENDDO
	545	DO k = 1, nz-1
	546	DO i = 0, nx
	547	tri(5,i,k) = tri(2,i,k) / tri(4,i,k-1)
	548	tri(4,i,k) = tri(1,i,k) - tri(3,i,k-1) * tri(5,i,k)
	549	ENDDO
	550	ENDDO
	551
	552	END SUBROUTINE split_1dd
	553
	554
	555	SUBROUTINE substi_1dd( ar, tri )
	556
	557	!------------------------------------------------------------------------------!
	558	! Substitution (Forward and Backward) (Thomas algorithm)
	559	!------------------------------------------------------------------------------!
	560
	561	IMPLICIT NONE
	562
	563	INTEGER :: i, k
	564
	565	REAL, DIMENSION(0:nx,nz) :: ar
	566	REAL, DIMENSION(0:nx,0:nz-1) :: ar1
	567	REAL, DIMENSION(5,0:nx,0:nz-1) :: tri
	568
	569	!
	570	!-- Forward substitution
	571	DO i = 0, nx
	572	ar1(i,0) = ar(i,1)
	573	ENDDO
	574	DO k = 1, nz-1
	575	DO i = 0, nx
	576	ar1(i,k) = ar(i,k+1) - tri(5,i,k) * ar1(i,k-1)
	577	ENDDO
	578	ENDDO
	579
	580	!
	581	!-- Backward substitution
	582	!-- Note, the add of 1.0E-20 in the denominator is due to avoid divisions
	583	!-- by zero appearing if the pressure bc is set to neumann at the top of
	584	!-- the model domain.
	585	DO i = 0, nx
	586	ar(i,nz) = ar1(i,nz-1) / ( tri(4,i,nz-1) + 1.0E-20 )
	587	ENDDO
	588	DO k = nz-2, 0, -1
	589	DO i = 0, nx
	590	ar(i,k+1) = ( ar1(i,k) - tri(3,i,k) * ar(i,k+2) ) &
	591	/ tri(4,i,k)
	592	ENDDO
	593	ENDDO
	594
	595	!
	596	!-- Indices i=0, j=0 correspond to horizontally averaged pressure.
	597	!-- The respective values of ar should be zero at all k-levels if
	598	!-- acceleration of horizontally averaged vertical velocity is zero.
	599	IF ( ibc_p_b == 1 .AND. ibc_p_t == 1 ) THEN
	600	IF ( j == 0 ) THEN
	601	DO k = 1, nz
	602	ar(0,k) = 0.0
	603	ENDDO
	604	ENDIF
	605	ENDIF
	606
	607	END SUBROUTINE substi_1dd
	608
	609	END SUBROUTINE tridia_1dd
	610
	611
	612	END MODULE tridia_solver

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |