1 | subroutine matvec(u, v, ipar) |
---|
2 | |
---|
3 | ! routine to apply the discrete differential operator A |
---|
4 | ! v <-- Au |
---|
5 | ! A is the compact representation of transformed Laplacian operator |
---|
6 | ! |
---|
7 | ! Interior points: |
---|
8 | ! V <- A[U] = 1/det_J*(Cxx U_x)_x + 1/y_eta*(Cyy U_y)_y + 1/det_J*(Czz U_z)_z |
---|
9 | ! |
---|
10 | ! Boundary points: |
---|
11 | ! V <- B[U] |
---|
12 | ! B = boundary operator |
---|
13 | ! for Dirichlet conditions B = I |
---|
14 | ! for Neumann conditions B = |d/dn| |
---|
15 | ! where n is outward facing unit normal |
---|
16 | !********************************************************** |
---|
17 | |
---|
18 | !grid _info contains Grid(0:nlevels-1), current_grid_level |
---|
19 | use grid_info |
---|
20 | use intermediate_variables |
---|
21 | use boundary_information |
---|
22 | use mpi_parameters |
---|
23 | use counters_flags_etc |
---|
24 | implicit none |
---|
25 | include 'mpif.h' |
---|
26 | |
---|
27 | real :: U(Grid(current_grid_level)%n1_it,Grid(current_grid_level)%n2_it,Grid(current_grid_level)%n3_it) |
---|
28 | real :: V(Grid(current_grid_level)%n1_it,Grid(current_grid_level)%n2_it,Grid(current_grid_level)%n3_it) |
---|
29 | |
---|
30 | !Local variables |
---|
31 | real, dimension(:,:,:), allocatable :: dUdzeta,dUdxi,dUdeta |
---|
32 | real, dimension(:,:,:,:), allocatable :: scratch |
---|
33 | |
---|
34 | integer :: ipar(13),ii,flag,npts |
---|
35 | integer :: status(MPI_STATUS_SIZE) |
---|
36 | real :: h |
---|
37 | real :: end_values(2) |
---|
38 | integer :: nx,ny,nz,locnz |
---|
39 | real :: dx,dy,dz |
---|
40 | real :: xx, mean_value |
---|
41 | |
---|
42 | V=0.0 |
---|
43 | ii = current_grid_level |
---|
44 | nx = Grid(ii)%n1 |
---|
45 | ny = Grid(ii)%n2 |
---|
46 | nz = Grid(ii)%n3 |
---|
47 | locnz = Grid(ii)%locn3 |
---|
48 | dx = Grid(ii)%dxi |
---|
49 | dy = Grid(ii)%deta |
---|
50 | dz = Grid(ii)%dzeta |
---|
51 | |
---|
52 | flag=1 ! don't specify deriv values at ends |
---|
53 | end_values=-999999. ! obsolete but required |
---|
54 | |
---|
55 | |
---|
56 | !*************************************************************************** |
---|
57 | ! Store approximation of A[U] in V: |
---|
58 | ! V <-- 1/det_J*(Cxx U_x)_x + 1/y_eta*(Cyy U_y)_y + 1/det_J*(Czz U_z)_z |
---|
59 | !*************************************************************************** |
---|
60 | |
---|
61 | |
---|
62 | !****************************************************************** |
---|
63 | ! Logic branches |
---|
64 | ! depending on xy periodicity |
---|
65 | !****************************************************************** |
---|
66 | |
---|
67 | if( boundary_type2(1,1)=='periodic' .and. boundary_type4(1,1)=='periodic' ) then |
---|
68 | |
---|
69 | allocate ( dUdzeta(1,1,2) ) |
---|
70 | |
---|
71 | !*************************************************************************** |
---|
72 | ! z part of the problem first, U & V are nz vectors |
---|
73 | !*************************************************************************** |
---|
74 | |
---|
75 | h=dz |
---|
76 | call compact_ddz(U,V,flag,end_values,1,1,nz) |
---|
77 | |
---|
78 | ! Save dU/dzeta |
---|
79 | if( myid==0 .and. BC_pressure5(1,1) == 'Neumann' ) then |
---|
80 | dUdzeta(1,1,1) = V(1,1,1) |
---|
81 | endif |
---|
82 | if( myid==nprocs-1 .and. BC_pressure6(1,1) == 'Neumann' ) then |
---|
83 | dUdzeta(1,1,2) = V(1,1,nz) |
---|
84 | endif |
---|
85 | |
---|
86 | !multiply U_z by Czz overwriting the result in held in V |
---|
87 | V(1,1,:) = Grid(ii)%Czz(1,:)*V(1,1,:) ! Czz is indep of x |
---|
88 | |
---|
89 | !differentiate (Czz U_z ) wrt z, overwrite results back into V |
---|
90 | call compact_ddz(V,V,flag,end_values,1,1,nz) |
---|
91 | |
---|
92 | !multiply result by 1/det_J |
---|
93 | V(1,1,:)=V(1,1,:)/Grid(ii)%det_J(1,:) |
---|
94 | |
---|
95 | ! We've aready done the z part of the operator, now add on the algebraic, horizontal part. |
---|
96 | ! We're working here in xy transform space so U is the current xy transformed P(z) estimate |
---|
97 | ! at a fixed kx,ky location (kx_val,ky_val; available via counters_flags_etc) |
---|
98 | |
---|
99 | V(1,1,:) = V(1,1,:) & |
---|
100 | -( (kx_val)**2 * (Grid(ii)%Cxx(1,:)/Grid(ii)%det_J(1,:)) & |
---|
101 | +(ky_val)**2 * (Grid(ii)%Cyy(1)/Grid(ii)%y_eta(1)) ) * U(1,1,:) |
---|
102 | |
---|
103 | ! Neumann conditions at top&bottom |
---|
104 | V(1,1,1)=dUdzeta(1,1,1) |
---|
105 | V(1,1,nz)=dUdzeta(1,1,2) |
---|
106 | |
---|
107 | deallocate (dUdzeta ) |
---|
108 | |
---|
109 | else ! NONPERIODIC LOGIC |
---|
110 | |
---|
111 | !*************************************************************************** |
---|
112 | ! z part of the problem first |
---|
113 | !*************************************************************************** |
---|
114 | allocate ( dUdzeta(nx,ny,2), dUdxi(ny,locnz,2), dUdeta(nx,locnz,2) ) |
---|
115 | |
---|
116 | h=dz |
---|
117 | call compact_ddz_mpi(U,V,flag,end_values,nx,ny,locnz) |
---|
118 | |
---|
119 | ! Save dU/dzeta |
---|
120 | if( myid==0 .and. BC_pressure5(1,1) == 'Neumann' ) then |
---|
121 | dUdzeta(:,:,1) = V(:,:,1) |
---|
122 | endif |
---|
123 | if( myid==nprocs-1 .and. BC_pressure6(1,1) == 'Neumann' ) then |
---|
124 | dUdzeta(:,:,2) = V(:,:,locnz) |
---|
125 | endif |
---|
126 | |
---|
127 | !multiply U_z by Czz overwriting the result in held in V |
---|
128 | do j=1,ny |
---|
129 | V(:,j,:) = Grid(ii)%Czz*V(:,j,:) |
---|
130 | enddo |
---|
131 | |
---|
132 | !differentiate (Czz U_z ) wrt z, overwrite results back into V |
---|
133 | call compact_ddz_mpi(V,V,flag,end_values,nx,ny,locnz) |
---|
134 | |
---|
135 | !multiply result by 1/det_J |
---|
136 | do j=1,ny |
---|
137 | V(:,j,:)=V(:,j,:)/Grid(ii)%det_J |
---|
138 | enddo |
---|
139 | |
---|
140 | !****************************************************************** |
---|
141 | ! do x part, compact differentiation |
---|
142 | !****************************************************************** |
---|
143 | |
---|
144 | !Take x derivative of U(:,:,1:locnz) and store result in scratch(:,:,:,1), |
---|
145 | h=dx |
---|
146 | call compact_ddx(U,scratch(1,1,1,1),flag,end_values,nx,ny,locnz) |
---|
147 | |
---|
148 | ! Save dU/dxi on FACES 1&2 |
---|
149 | dUdxi(:,:,1) = scratch(1,:,:,1) |
---|
150 | dUdxi(:,:,2) = scratch(nx,:,:,1) |
---|
151 | |
---|
152 | !multiply U_x by Cxx, overwriting the result in scratch(:,:,:,1) |
---|
153 | do j=1,ny |
---|
154 | scratch(:,j,:,1) = Grid(ii)%Cxx*scratch(:,j,:,1) |
---|
155 | enddo |
---|
156 | |
---|
157 | !differentiate (Cxx U_x) wrt x, writing results into scratch(:,:,:,2) |
---|
158 | call compact_ddx(scratch(1,1,1,1),scratch(1,1,1,2),flag,end_values,nx,ny,locnz) |
---|
159 | |
---|
160 | !divide (Cxx U_x)_x by det_J and add the result to V |
---|
161 | !V then contains 1/det_J*{ (Cxx U_x)_x + (Czz U_z)_z } |
---|
162 | do j=1,ny |
---|
163 | V(:,j,:) = V(:,j,:) + scratch(:,j,:,2)/Grid(ii)%det_J |
---|
164 | enddo |
---|
165 | |
---|
166 | !****************************************************************** |
---|
167 | ! do y part, compact differentiation |
---|
168 | !****************************************************************** |
---|
169 | |
---|
170 | !Take y derivative of U(:,:,1:locnz) and store result in scratch(:,:,:,1). |
---|
171 | h=dy |
---|
172 | call compact_ddy(U,scratch(1,1,1,1),flag,end_values,nx,ny,locnz) |
---|
173 | |
---|
174 | ! Save dU/deta on FACES 3&4 |
---|
175 | dUdeta(:,:,1) = scratch(:,1,:,1) |
---|
176 | dUdeta(:,:,2) = scratch(:,ny,:,1) |
---|
177 | |
---|
178 | !multiply U_y by Cyy, overwriting the result in scratch |
---|
179 | do j=1,ny |
---|
180 | scratch(:,j,:,1) = Grid(ii)%Cyy(j)*scratch(:,j,:,1) |
---|
181 | enddo |
---|
182 | |
---|
183 | !differentiate (Cyy U_y) wrt y, write results back into scratch(:,:,:,2) |
---|
184 | call compact_ddy(scratch(1,1,1,1),scratch(1,1,1,2),flag,end_values,nx,ny,locnz) |
---|
185 | |
---|
186 | !add 1/y_eta (Cyy U_y)_y to V. Note Cyy=1/y_eta. V then contains: |
---|
187 | !1/det_J*(Cxx U_x)_x + 1/y_eta*(Cyy U_y)_y + 1/det_J*(Czz U_z)_z |
---|
188 | do j=1,ny |
---|
189 | V(:,j,:) = V(:,j,:) + Grid(ii)%Cyy(j)*scratch(:,j,:,2) |
---|
190 | enddo |
---|
191 | |
---|
192 | !****************************************************************** |
---|
193 | ! overwrite boundary locations w/ boundary operators |
---|
194 | !****************************************************************** |
---|
195 | |
---|
196 | ! USE IDENTITY OPERATOR FOR DIRICHLET BCs:************************* |
---|
197 | if( BC_pressure1(1,1) == 'Dirichlet' ) then |
---|
198 | V(1,:,:)=U(1,:,:) |
---|
199 | endif |
---|
200 | if( BC_pressure2(1,1) == 'Dirichlet' ) then |
---|
201 | V(nx,:,:)=U(nx,:,:) |
---|
202 | endif |
---|
203 | |
---|
204 | if( BC_pressure3(1,1) == 'Dirichlet' ) then |
---|
205 | V(:,1,:)=U(:,1,:) |
---|
206 | endif |
---|
207 | if( BC_pressure4(1,1) == 'Dirichlet' ) then |
---|
208 | V(:,ny,:)=U(:,ny,:) |
---|
209 | endif |
---|
210 | |
---|
211 | if( myid==0 .and. BC_pressure5(1,1) == 'Dirichlet' ) then |
---|
212 | V(:,:,1)=U(:,:,1) |
---|
213 | endif |
---|
214 | if( myid==nprocs-1 .and. BC_pressure6(1,1) == 'Dirichlet' ) then |
---|
215 | V(:,:,locnz)=U(:,:,locnz) |
---|
216 | endif |
---|
217 | |
---|
218 | ! USE NORMAL DERIV. OPERATOR FOR NEUMANN BCs:************************* |
---|
219 | if( BC_pressure1(1,1) == 'Neumann' ) then |
---|
220 | V(1,:,:)=dUdxi(:,:,1) |
---|
221 | endif |
---|
222 | if( BC_pressure2(1,1) == 'Neumann' ) then |
---|
223 | V(nx,:,:)=dUdxi(:,:,2) |
---|
224 | endif |
---|
225 | |
---|
226 | if( BC_pressure3(1,1) == 'Neumann' ) then |
---|
227 | V(:,1,:)=dUdeta(:,:,1) |
---|
228 | endif |
---|
229 | if( BC_pressure4(1,1) == 'Neumann' ) then |
---|
230 | V(:,ny,:)=dUdeta(:,:,2) |
---|
231 | endif |
---|
232 | |
---|
233 | if( myid==0 .and. BC_pressure5(1,1) == 'Neumann' ) then |
---|
234 | V(:,:,1)=dUdzeta(:,:,1) |
---|
235 | endif |
---|
236 | if( myid==nprocs-1 .and. BC_pressure6(1,1) == 'Neumann' ) then |
---|
237 | V(:,:,locnz)=dUdzeta(:,:,2) |
---|
238 | endif |
---|
239 | |
---|
240 | deallocate (dUdzeta, dUdxi,dUdeta) |
---|
241 | endif ! END XY PERIODIC IF BLOCK |
---|
242 | |
---|
243 | end subroutine matvec |
---|
244 | |
---|
245 | |
---|
246 | |
---|
247 | subroutine mg_solve(b,x,ipar) |
---|
248 | |
---|
249 | ! Routine to apply M ~ A_inverse to the rhs |
---|
250 | ! vector b, producing an approximate solution x. |
---|
251 | ! M is defined to be n_mg_cycles full multigrid |
---|
252 | ! cycles. The FMG cycles start and end at the |
---|
253 | ! finest (h=0) grid level. |
---|
254 | ! |
---|
255 | |
---|
256 | use grid_info |
---|
257 | use counters_flags_etc, only: n_mg_cycles,istep,istart,tol |
---|
258 | use mpi_parameters, only: myid |
---|
259 | implicit none |
---|
260 | |
---|
261 | integer :: i,N,ipar(13) |
---|
262 | real :: b(Grid(0)%n1,Grid(0)%n2,Grid(0)%locn3) |
---|
263 | real :: x(Grid(0)%n1,Grid(0)%n2,Grid(0)%locn3) |
---|
264 | character(len=3) :: mg_cycle_type |
---|
265 | |
---|
266 | if( istep<istart ) then |
---|
267 | mg_cycle_type='FMG' |
---|
268 | else |
---|
269 | mg_cycle_type='V' |
---|
270 | endif |
---|
271 | N=n_mg_cycles |
---|
272 | |
---|
273 | ! (1) store the input vector b in the "f" location of the fine grid: |
---|
274 | Grid(0)%f = b |
---|
275 | |
---|
276 | ! (2) MG cycles: |
---|
277 | if( mg_cycle_type=='FMG') then |
---|
278 | call FMG(0) ! use Grid(0)%f as input and updates Grid(0)%v |
---|
279 | elseif(mg_cycle_type=='V') then |
---|
280 | Grid(0)%v=0.0 |
---|
281 | do i=1,N |
---|
282 | call V_cycle(0) |
---|
283 | if(myid==0)write(0,*) i,Grid(0)%spar(2) |
---|
284 | if( Grid(0)%spar(2) < tol ) goto 999 ! convergence |
---|
285 | enddo |
---|
286 | else |
---|
287 | write(0,*) 'mg_cycle_type not properly defined' |
---|
288 | STOP |
---|
289 | endif |
---|
290 | |
---|
291 | ! (3) extract the solution from Grid(0), store in output vector x: |
---|
292 | 999 continue |
---|
293 | x = Grid(0)%v |
---|
294 | |
---|
295 | current_grid_level=0 ! return global variable to fine grid value |
---|
296 | end subroutine mg_solve |
---|
297 | |
---|
298 | |
---|
299 | subroutine diagl(U,V,ipar) |
---|
300 | ! |
---|
301 | ! Routine to apply M ~ (1/diag) to the rhs |
---|
302 | ! vector b, producing an approximate solution x. |
---|
303 | ! |
---|
304 | |
---|
305 | use grid_info |
---|
306 | use boundary_information |
---|
307 | use counters_flags_etc |
---|
308 | use mpi_parameters |
---|
309 | implicit none |
---|
310 | |
---|
311 | integer :: h,ipar(13) |
---|
312 | real :: diag_FD ! approximate diagonal element in a FD discretization at level h |
---|
313 | real :: pi |
---|
314 | real :: U(Grid(current_grid_level)%n1,Grid(current_grid_level)%n2,Grid(current_grid_level)%locn3) |
---|
315 | real :: V(Grid(current_grid_level)%n1,Grid(current_grid_level)%n2,Grid(current_grid_level)%locn3) |
---|
316 | real, allocatable :: kx(:),ky(:) |
---|
317 | integer :: N |
---|
318 | |
---|
319 | |
---|
320 | h=current_grid_level |
---|
321 | |
---|
322 | if( boundary_type2(1,1)=='periodic' .and. boundary_type2(1,1)=='periodic' ) then |
---|
323 | |
---|
324 | ! create the x & y wavenumber arrays, align with the data arrays |
---|
325 | N=Grid(h)%n1-1 ! n1 includes both boundary points, N only 1 |
---|
326 | pi=4.*atan(1.) |
---|
327 | allocate(kx(N+1)) |
---|
328 | do i=1,N/2+1 |
---|
329 | kx(i)=2.*pi*(i-1.) ! FFTW real to half-complex storage |
---|
330 | enddo |
---|
331 | do i=2,N/2 |
---|
332 | kx(N-i+2)=kx(i) ! FFTW real to half-complex storage |
---|
333 | enddo |
---|
334 | kx(N+1)=0. ! extra location, convenient for alignment of arrays |
---|
335 | |
---|
336 | |
---|
337 | N=Grid(h)%n2-1 ! n2 includes both boundary points, N only 1 |
---|
338 | allocate(ky(N+1)) |
---|
339 | do j=1,N/2+1 |
---|
340 | ky(j)=2.*pi*(j-1.) ! FFTW real to half-complex storage |
---|
341 | enddo |
---|
342 | do j=2,N/2 |
---|
343 | ky(N-j+2)=ky(j) ! FFTW real to half-complex storage |
---|
344 | enddo |
---|
345 | ky(N+1)=0. ! extra location, convenient for alignment of arrays |
---|
346 | |
---|
347 | ! precondition using the diagonal elements of an approximated operator |
---|
348 | do i=1,Grid(h)%n1 |
---|
349 | do j=1,Grid(h)%n2 |
---|
350 | do k=1,Grid(h)%locn3 |
---|
351 | diag_FD = -2.*( Grid(h)%Czz(i,k)/(Grid(h)%dzeta)**2 ) & ! diagonal term for FD approx of z operator |
---|
352 | - (kx(i)**2) * (Grid(h)%Cxx(i,k)/(Grid(h)%det_J(i,k))) & ! true diag term in F space |
---|
353 | - (ky(j)**2) * (Grid(h)%Cyy(j)/(Grid(h)%y_eta(j))) ! true diag term in F space |
---|
354 | |
---|
355 | if( j==1 .and. k==1 ) then |
---|
356 | write(0,*)'****** diag: ',i, -2.*( Grid(h)%Czz(i,k)/(Grid(h)%dzeta)**2 ), - (kx(i)**2) * (Grid(h)%Cxx(i,k)/(Grid(h)%det_J(i,k))) |
---|
357 | endif |
---|
358 | V(i,j,k) = U(i,j,k)/diag_FD ! interior points |
---|
359 | enddo |
---|
360 | enddo |
---|
361 | enddo |
---|
362 | STOP |
---|
363 | |
---|
364 | else |
---|
365 | |
---|
366 | ! precondition using the diagonal elements of an approximated operator |
---|
367 | do i=1,Grid(h)%n1 |
---|
368 | do j=1,Grid(h)%n2 |
---|
369 | do k=1,Grid(h)%locn3 |
---|
370 | diag_FD = -2.*( Grid(h)%Cxx(i,k)/(Grid(h)%dxi)**2 & |
---|
371 | + Grid(h)%Cyy(j )/(Grid(h)%deta)**2 & |
---|
372 | + Grid(h)%Czz(i,k)/(Grid(h)%dzeta)**2 ) |
---|
373 | V(i,j,k) = U(i,j,k)/diag_FD ! interior points |
---|
374 | enddo |
---|
375 | enddo |
---|
376 | enddo |
---|
377 | |
---|
378 | endif |
---|
379 | |
---|
380 | ! Redo Neumann boundaries: i.e. overwrite locations where the boundary operator differs from the interior operator |
---|
381 | |
---|
382 | if( BC_pressure1(1,1) == 'Neumann' ) then |
---|
383 | do k=1,Grid(h)%locn3 |
---|
384 | do j=1,Grid(h)%n2 |
---|
385 | diag_FD = 1./Grid(h)%dxi |
---|
386 | V(1,j,k)= U(1,j,k)/diag_FD |
---|
387 | enddo |
---|
388 | enddo |
---|
389 | endif |
---|
390 | |
---|
391 | if( BC_pressure2(1,1) == 'Neumann' ) then |
---|
392 | do k=1,Grid(h)%locn3 |
---|
393 | do j=1,Grid(h)%n2 |
---|
394 | diag_FD = 1./Grid(h)%dxi |
---|
395 | V(Grid(h)%n1,j,k)= U(Grid(h)%n1,j,k)/diag_FD |
---|
396 | enddo |
---|
397 | enddo |
---|
398 | endif |
---|
399 | |
---|
400 | if( BC_pressure3(1,1) == 'Neumann' ) then |
---|
401 | do k=1,Grid(h)%locn3 |
---|
402 | do i=1,Grid(h)%n1 |
---|
403 | diag_FD = 1./Grid(h)%deta |
---|
404 | V(i,1,k)= U(i,1,k)/diag_FD |
---|
405 | enddo |
---|
406 | enddo |
---|
407 | endif |
---|
408 | |
---|
409 | if( BC_pressure4(1,1) == 'Neumann' ) then |
---|
410 | do k=1,Grid(h)%locn3 |
---|
411 | do i=1,Grid(h)%n1 |
---|
412 | diag_FD = 1./Grid(h)%deta |
---|
413 | V(i,Grid(h)%n2,k)= U(i,Grid(h)%n2,k)/diag_FD |
---|
414 | enddo |
---|
415 | enddo |
---|
416 | endif |
---|
417 | |
---|
418 | if( BC_pressure5(1,1) == 'Neumann' .and. myid==0 ) then |
---|
419 | do i=1,Grid(h)%n1 |
---|
420 | do j=1,Grid(h)%n2 |
---|
421 | diag_FD = 1./Grid(h)%dzeta |
---|
422 | V(i,j,1)= U(i,j,1)/diag_FD |
---|
423 | enddo |
---|
424 | enddo |
---|
425 | endif |
---|
426 | |
---|
427 | if( BC_pressure6(1,1) == 'Neumann' .and. myid==nprocs-1) then |
---|
428 | do i=1,Grid(h)%n1 |
---|
429 | do j=1,Grid(h)%n2 |
---|
430 | diag_FD = 1./Grid(h)%dzeta |
---|
431 | V(i,j,Grid(h)%locn3)= U(i,j,Grid(h)%locn3)/diag_FD |
---|
432 | enddo |
---|
433 | enddo |
---|
434 | endif |
---|
435 | |
---|
436 | if( boundary_type2(1,1)=='periodic' ) then |
---|
437 | do k=1,Grid(h)%locn3 |
---|
438 | do j=1,Grid(h)%n2 |
---|
439 | diag_FD = 1. ! Identity operator |
---|
440 | V(Grid(h)%n1,j,k)= U(Grid(h)%n1,j,k)/diag_FD ! only "far" endpts have I[P]=0 |
---|
441 | enddo |
---|
442 | enddo |
---|
443 | endif |
---|
444 | |
---|
445 | if( boundary_type4(1,1)=='periodic' ) then |
---|
446 | do k=1,Grid(h)%locn3 |
---|
447 | do i=1,Grid(h)%n1 |
---|
448 | diag_FD = 1 |
---|
449 | V(i,Grid(h)%n2,k)= U(i,Grid(h)%n2,k)/diag_FD |
---|
450 | enddo |
---|
451 | enddo |
---|
452 | endif |
---|
453 | |
---|
454 | return |
---|
455 | end |
---|