Home

Context Navigation

sor.f90 @ 707

Last change on this file since 707 was 707, checked in by raasch, 13 years ago

New:
---

In case of multigrid method, on coarse grid levels, gathered data are
identically processed on all PEs (before, on PE0 only), so that the subsequent
scattering of data is not neccessary any more. (modules, init_pegrid, poismg)

Changed:

Calculation of weighted average of p is now handled in the same way
regardless of the number of ghost layers (advection scheme). (pres)

multigrid and sor method are using p_loc for iterative
advancements of pressure. p_sub removed. (init_3d_model, modules, poismg, pres, sor)

bc_lr and bc_ns replaced by bc_lr_dirrad, bc_lr_raddir, bc_ns_dirrad, bc_ns_raddir
for speed optimization. (calc_spectra, check_parameters, exchange_horiz,
exchange_horiz_2d, header, init_3d_model, init_grid, init_pegrid, modules,
poismg, pres, sor, time_integration, timestep)

grid_level directly used as index for MPI data type arrays. (exchange_horiz,
poismg)

initial assignments of zero to array p for iterative solvers only (init_3d_model)

Errors:

localsum calculation modified for proper OpenMP reduction. (pres)

Bugfix: bottom (nzb) and top (nzt+1) boundary conditions set in routines
resid and restrict. They were missed before, which may have led to
unpredictable results. (poismg)

Property svn:keywords set to Id

File size: 6.0 KB

Line
1	SUBROUTINE sor( d, ddzu, ddzw, p )
2
3	!------------------------------------------------------------------------------!
4	! Current revisions:
5	! -----------------
6	! bc_lr/ns replaced by bc_lr/ns_cyc
7	!
8	! Former revisions:
9	! -----------------
10	! $Id: sor.f90 707 2011-03-29 11:39:40Z raasch $
11	!
12	! 667 2010-12-23 12:06:00Z suehring/gryschka
13	! nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng.
14	! Call of exchange_horiz are modified.
15	! bug removed in declaration of ddzw(), nz replaced by nzt+1
16	!
17	! 75 2007-03-22 09:54:05Z raasch
18	! 2nd+3rd argument removed from exchange horiz
19	!
20	! RCS Log replace by Id keyword, revision history cleaned up
21	!
22	! Revision 1.9 2005/03/26 21:02:23 raasch
23	! Implementation of non-cyclic (Neumann) horizontal boundary conditions,
24	! dx2,dy2 replaced by ddx2,ddy2
25	!
26	! Revision 1.1 1997/08/11 06:25:56 raasch
27	! Initial revision
28	!
29	!
30	! Description:
31	! ------------
32	! Solve the Poisson-equation with the SOR-Red/Black-scheme.
33	!------------------------------------------------------------------------------!
34
35	USE grid_variables
36	USE indices
37	USE pegrid
38	USE control_parameters
39
40	IMPLICIT NONE
41
42	INTEGER :: i, j, k, n, nxl1, nxl2, nys1, nys2
43	REAL :: ddzu(1:nz+1), ddzw(1:nzt+1)
44	REAL :: d(nzb+1:nzt,nys:nyn,nxl:nxr), &
45	p(nzb:nzt+1,nysg:nyng,nxlg:nxrg)
46	REAL, DIMENSION(:), ALLOCATABLE :: f1, f2, f3
47
48	ALLOCATE( f1(1:nz), f2(1:nz), f3(1:nz) )
49
50	!
51	!-- Compute pre-factors.
52	DO k = 1, nz
53	f2(k) = ddzu(k+1) * ddzw(k)
54	f3(k) = ddzu(k) * ddzw(k)
55	f1(k) = 2.0 * ( ddx2 + ddy2 ) + f2(k) + f3(k)
56	ENDDO
57
58	!
59	!-- Limits for RED- and BLACK-part.
60	IF ( MOD( nxl , 2 ) == 0 ) THEN
61	nxl1 = nxl
62	nxl2 = nxl + 1
63	ELSE
64	nxl1 = nxl + 1
65	nxl2 = nxl
66	ENDIF
67	IF ( MOD( nys , 2 ) == 0 ) THEN
68	nys1 = nys
69	nys2 = nys + 1
70	ELSE
71	nys1 = nys + 1
72	nys2 = nys
73	ENDIF
74
75	DO n = 1, n_sor
76
77	!
78	!-- RED-part
79	DO i = nxl1, nxr, 2
80	DO j = nys2, nyn, 2
81	DO k = nzb+1, nzt
82	p(k,j,i) = p(k,j,i) + omega_sor / f1(k) * ( &
83	ddx2 * ( p(k,j,i+1) + p(k,j,i-1) ) + &
84	ddy2 * ( p(k,j+1,i) + p(k,j-1,i) ) + &
85	f2(k) * p(k+1,j,i) + &
86	f3(k) * p(k-1,j,i) - &
87	d(k,j,i) - &
88	f1(k) * p(k,j,i) )
89	ENDDO
90	ENDDO
91	ENDDO
92
93	DO i = nxl2, nxr, 2
94	DO j = nys1, nyn, 2
95	DO k = nzb+1, nzt
96	p(k,j,i) = p(k,j,i) + omega_sor / f1(k) * ( &
97	ddx2 * ( p(k,j,i+1) + p(k,j,i-1) ) + &
98	ddy2 * ( p(k,j+1,i) + p(k,j-1,i) ) + &
99	f2(k) * p(k+1,j,i) + &
100	f3(k) * p(k-1,j,i) - &
101	d(k,j,i) - &
102	f1(k) * p(k,j,i) )
103	ENDDO
104	ENDDO
105	ENDDO
106
107	!
108	!-- Exchange of boundary values for p.
109	CALL exchange_horiz( p, nbgp )
110
111	!
112	!-- Horizontal (Neumann) boundary conditions in case of non-cyclic boundaries
113	IF ( .NOT. bc_lr_cyc ) THEN
114	IF ( inflow_l .OR. outflow_l ) p(:,:,nxl-1) = p(:,:,nxl)
115	IF ( inflow_r .OR. outflow_r ) p(:,:,nxr+1) = p(:,:,nxr)
116	ENDIF
117	IF ( .NOT. bc_ns_cyc ) THEN
118	IF ( inflow_n .OR. outflow_n ) p(:,nyn+1,:) = p(:,nyn,:)
119	IF ( inflow_s .OR. outflow_s ) p(:,nys-1,:) = p(:,nys,:)
120	ENDIF
121
122	!
123	!-- BLACK-part
124	DO i = nxl1, nxr, 2
125	DO j = nys1, nyn, 2
126	DO k = nzb+1, nzt
127	p(k,j,i) = p(k,j,i) + omega_sor / f1(k) * ( &
128	ddx2 * ( p(k,j,i+1) + p(k,j,i-1) ) + &
129	ddy2 * ( p(k,j+1,i) + p(k,j-1,i) ) + &
130	f2(k) * p(k+1,j,i) + &
131	f3(k) * p(k-1,j,i) - &
132	d(k,j,i) - &
133	f1(k) * p(k,j,i) )
134	ENDDO
135	ENDDO
136	ENDDO
137
138	DO i = nxl2, nxr, 2
139	DO j = nys2, nyn, 2
140	DO k = nzb+1, nzt
141	p(k,j,i) = p(k,j,i) + omega_sor / f1(k) * ( &
142	ddx2 * ( p(k,j,i+1) + p(k,j,i-1) ) + &
143	ddy2 * ( p(k,j+1,i) + p(k,j-1,i) ) + &
144	f2(k) * p(k+1,j,i) + &
145	f3(k) * p(k-1,j,i) - &
146	d(k,j,i) - &
147	f1(k) * p(k,j,i) )
148	ENDDO
149	ENDDO
150	ENDDO
151
152	!
153	!-- Exchange of boundary values for p.
154	CALL exchange_horiz( p, nbgp )
155
156	!
157	!-- Boundary conditions top/bottom.
158	!-- Bottom boundary
159	IF ( ibc_p_b == 1 ) THEN ! Neumann
160	p(nzb,:,:) = p(nzb+1,:,:)
161	ELSE ! Dirichlet
162	p(nzb,:,:) = 0.0
163	ENDIF
164
165	!
166	!-- Top boundary
167	IF ( ibc_p_t == 1 ) THEN ! Neumann
168	p(nzt+1,:,:) = p(nzt,:,:)
169	ELSE ! Dirichlet
170	p(nzt+1,:,:) = 0.0
171	ENDIF
172
173	!
174	!-- Horizontal (Neumann) boundary conditions in case of non-cyclic boundaries
175	IF ( .NOT. bc_lr_cyc ) THEN
176	IF ( inflow_l .OR. outflow_l ) p(:,:,nxl-1) = p(:,:,nxl)
177	IF ( inflow_r .OR. outflow_r ) p(:,:,nxr+1) = p(:,:,nxr)
178	ENDIF
179	IF ( .NOT. bc_ns_cyc ) THEN
180	IF ( inflow_n .OR. outflow_n ) p(:,nyn+1,:) = p(:,nyn,:)
181	IF ( inflow_s .OR. outflow_s ) p(:,nys-1,:) = p(:,nys,:)
182	ENDIF
183
184
185	ENDDO
186
187	DEALLOCATE( f1, f2, f3 )
188
189	END SUBROUTINE sor

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

Context Navigation

source: palm/trunk/SOURCE/sor.f90 @ 707

Download in other formats: