20 procedure(cg),
pointer:: update => null()
25 procedure,
nopass:: laplace
27 real,
save:: floor=1e-2, tolerance=1e-4, fourpig=1.0
28 integer,
save:: verbose=0
29 integer,
save:: max_iter=500
30 logical,
save:: precondition=.true.
31 real(8),
parameter:: pi = asin(1.0) * 2.0
32 logical,
save:: chebyshev = .true.
39 SUBROUTINE init (self)
45 character(len=12),
save :: solver=
'cg' 46 logical,
save :: first_time = .true.
47 namelist /poisson_params/ tolerance, floor, max_iter, fourpig, precondition, &
48 verbose, solver, chebyshev
53 fourpig = scaling%grav
55 read (io%input, poisson_params, iostat=iostat)
56 if (io%master)
write (io%output, poisson_params)
59 select case(trim(solver))
66 write(*,*)
'Poisson solver '//trim(solver)//
' is unknown' 67 write(*,*)
'Valid solvers are: cg, sor' 81 SUBROUTINE cg (self, patch, phi1, d, d0, detailed_timer)
84 real(kind=KindScalarVar),
dimension(:,:,:),
pointer:: phi1
85 real(kind=KindScalarVar),
dimension(:,:,:),
pointer,
intent(in):: d
87 logical :: detailed_timer
90 real :: error, error0, usppt
91 real(8) :: pap, ru, ru2, alpha, beta, ds(3), wc
92 integer :: iter, m(3), ione, jone, kone
93 integer,
save :: itimer=0
94 real(8),
dimension(:,:,:),
allocatable:: source, phi, r, p, s, u
96 subroutine laplace_interface(ds, phi, out)
97 real(8),
intent(in) :: ds(3)
98 real(8),
dimension(:,:,:),
intent(in) :: phi
99 real(8),
dimension(:,:,:),
intent(out):: out
100 end subroutine laplace_interface
102 procedure(laplace_interface),
pointer:: laplace
104 if (detailed_timer)
call trace%begin(
'poisson_mod%cg', itimer=itimer)
108 if (
size(phi1,1) <= 1) ione = 0
109 if (
size(phi1,2) <= 1) jone = 0
110 if (
size(phi1,3) <= 1) kone = 0
112 if (patch%mesh_type == mesh_types%Cartesian)
then 113 laplace => laplace_cartesian
115 laplace => laplace_curvilinear
121 allocate (r(m(1),m(2),m(3)), p(m(1),m(2),m(3)))
122 allocate (s(m(1),m(2),m(3)), u(m(1),m(2),m(3)))
123 allocate (source(m(1),m(2),m(3)), phi(m(1),m(2),m(3)))
130 source = fourpig*(d-d0)
131 call zero_boundaries (source)
136 call laplace (ds, phi, r)
138 error0 = patch%fmaxval(abs(r/(d+floor)))
139 if (verbose > 1) print
"(i6,2x,'cg: iter = ',i4.1,3x,'error =',1p,g12.4)", &
148 call laplace (ds, p, s)
149 alpha = ru/patch%fsum(s*p)
152 error = patch%fmaxval(abs(r/(d+floor)))/fourpig
154 ru2 = patch%fsum(r*u)
155 if (error < tolerance)
exit 159 if (verbose > 1) print
"(i6,2x,'cg: iter = ',i4.1,3x,'error =',1p,g12.4)", &
160 patch%id, iter, error
163 if (iter==max_iter+1)
then 164 print *,
'poisson_t%cg: NO CONVERGENCE! id, niter, error0, error =', &
165 patch%id, max_iter, error0, error
166 else if (verbose>0 .and. io%master)
then 167 usppt = 1e6*(wallclock()-wc)/product(shape(phi))
168 print
'(a,i7,i5,3x,2e12.3,f8.2)', &
169 'poisson_t%cg: id, niter, initial/final error, us/pt =', &
170 patch%id, iter, error0, error, usppt
174 deallocate (phi, source, r, p, s, u)
175 if (detailed_timer)
call trace%end (itimer)
180 subroutine zero_boundaries (f)
182 class(
mesh_t),
pointer:: m(:)
184 if (patch%is_periodic())
return 186 if (
size(f,1) > 1)
then 187 f(m(1)%lb:m(1)%lo,:,:) = 0d0
188 f(m(1)%uo:m(1)%ub,:,:) = 0d0
190 if (
size(f,2) > 1)
then 191 f(:,m(2)%lb:m(2)%lo,:) = 0d0
192 f(:,m(2)%uo:m(2)%ub,:) = 0d0
194 if (
size(f,3) > 1)
then 195 f(:,:,m(3)%lb:m(3)%lo) = 0d0
196 f(:,:,m(3)%uo:m(3)%ub) = 0d0
198 end subroutine zero_boundaries
202 subroutine laplace_cartesian (ds, phi, out)
203 real(8),
intent(in) :: ds(3)
204 real(8),
dimension(:,:,:),
intent(in) :: phi
205 real(8),
dimension(:,:,:),
intent(out):: out
207 integer :: i, j, k, im1, jm1, km1, ip1, jp1, kp1, li(3), n(3)
219 if (patch%is_periodic())
then 220 do k=patch%mesh(3)%lb,patch%mesh(3)%ub
221 kp1 = modulo(k-li(3)+1,n(3))+li(3)
222 km1 = modulo(k-li(3)-1,n(3))+li(3)
223 do j=patch%mesh(2)%lb,patch%mesh(2)%ub
224 jp1 = modulo(j-li(2)+1,n(2))+li(2)
225 jm1 = modulo(j-li(2)-1,n(2))+li(2)
226 do i=patch%mesh(1)%lb,patch%mesh(1)%lo
227 ip1 = modulo(i-li(1)+1,n(1))+li(1)
228 im1 = modulo(i-li(1)-1,n(1))+li(1)
229 out(i,j,k) = c(1)*(phi(ip1,j,k)+phi(im1,j,k)-2d0*phi(i,j,k)) &
230 + c(2)*(phi(i,jp1,k)+phi(i,jm1,k)-2d0*phi(i,j,k)) &
231 + c(3)*(phi(i,j,kp1)+phi(i,j,km1)-2d0*phi(i,j,k))
233 do i=patch%mesh(1)%li,patch%mesh(1)%ui
234 out(i,j,k) = c(1)*(phi(i+ione,j,k)+phi(i-ione,j,k)-2d0*phi(i,j,k)) &
235 + c(2)*(phi(i ,jp1,k)+phi(i ,jm1,k)-2d0*phi(i,j,k)) &
236 + c(3)*(phi(i ,j,kp1)+phi(i ,j,km1)-2d0*phi(i,j,k))
238 do i=patch%mesh(1)%uo,patch%mesh(1)%ub
239 ip1 = modulo(i-li(1)+1,n(1))+li(1)
240 im1 = modulo(i-li(1)-1,n(1))+li(1)
241 out(i,j,k) = c(1)*(phi(ip1,j,k)+phi(im1,j,k)-2d0*phi(i,j,k)) &
242 + c(2)*(phi(i,jp1,k)+phi(i,jm1,k)-2d0*phi(i,j,k)) &
243 + c(3)*(phi(i,j,kp1)+phi(i,j,km1)-2d0*phi(i,j,k))
251 call zero_boundaries (out)
252 do k=patch%mesh(3)%li,patch%mesh(3)%ui
253 do j=patch%mesh(2)%li,patch%mesh(2)%ui
254 do i=patch%mesh(1)%li,patch%mesh(1)%ui
255 out(i,j,k) = c(1)*(phi(i+ione,j,k)+phi(i-ione,j,k)-2d0*phi(i,j,k)) &
256 + c(2)*(phi(i,j+jone,k)+phi(i,j-jone,k)-2d0*phi(i,j,k)) &
257 + c(3)*(phi(i,j,k+kone)+phi(i,j,k-kone)-2d0*phi(i,j,k))
262 end subroutine laplace_cartesian
266 subroutine laplace_curvilinear (ds, phi, out)
267 real(8),
intent(in) :: ds(3)
268 real(8),
dimension(:,:,:),
intent(in) :: phi
269 real(8),
dimension(:,:,:),
intent(out):: out
271 class(
mesh_t),
pointer:: m1, m2
273 real(8) :: c(3), h2ci, h31ci, h32ci
276 where (patch%mesh%nc == 1)
286 if (patch%is_periodic())
then 287 call mpi%abort(
'NOT IMPLEMENTED ERROR: This curvilinear patch is periodic in all three directions. Really??')
292 call zero_boundaries (out)
293 do k=patch%mesh(3)%li,patch%mesh(3)%ui
294 do j=patch%mesh(2)%li,patch%mesh(2)%ui
295 h32ci = 1.0d0 / m2%h32c(j)
296 do i=patch%mesh(1)%li,patch%mesh(1)%ui
297 h2ci = 1.0d0 / m1%h2c(i)
298 h31ci = 1.0d0 / m1%h31c(i)
299 out(i,j,k) = c(1) * h2ci * h31ci &
300 * (m1%h2f(i+ione) * m1%h31f(i+ione) * (phi(i+ione,j,k) - phi(i,j,k)) &
301 - m1%h2f(i ) * m1%h31f(i ) * (phi(i,j,k) - phi(i-ione,j,k))) &
302 + c(2) * h2ci**2 * h32ci &
303 * (m2%h32f(j+jone) * (phi(i,j+jone,k) - phi(i,j,k)) &
304 - m2%h32f(j ) * (phi(i,j,k) - phi(i,j-jone,k))) &
305 + c(3) * h31ci**2 * h32ci**2 &
306 * (phi(i,j,k+kone) + phi(i,j,k-kone) - 2.0d0 * phi(i,j,k))
311 end subroutine laplace_curvilinear
315 subroutine precond (res, out)
316 real(8),
dimension(:,:,:) :: res, out
318 integer :: i, j, k, im1, jm1, km1, ip1, jp1, kp1, li(3), n(3), &
322 if (precondition)
then 324 ndim = sum(merge(1,0,n > 1))
326 a = 1d0 + 0.25d0 / ndim
334 if (patch%is_periodic())
then 335 do k=patch%mesh(3)%lb,patch%mesh(3)%ub
336 kp1 = modulo(k-li(3)+1,n(3))+li(3)
337 km1 = modulo(k-li(3)-1,n(3))+li(3)
338 do j=patch%mesh(2)%lb,patch%mesh(2)%ub
339 jp1 = modulo(j-li(2)+1,n(2))+li(2)
340 jm1 = modulo(j-li(2)-1,n(2))+li(2)
341 do i=patch%mesh(1)%lb,patch%mesh(1)%ub
342 ip1 = modulo(i-li(1)+1,n(1))+li(1)
343 im1 = modulo(i-li(1)-1,n(1))+li(1)
344 out(i,j,k) = b(1)*(res(ip1,j,k)+res(im1,j,k)) &
345 + b(2)*(res(i,jp1,k)+res(i,jm1,k)) &
346 + b(3)*(res(i,j,kp1)+res(i,j,km1)) &
355 do k=patch%li(3),patch%ui(3)
356 do j=patch%li(2),patch%ui(2)
357 do i=patch%li(1),patch%ui(1)
358 out(i,j,k) = b(1)*(res(i+ione,j,k)+res(i-ione,j,k)) &
359 + b(2)*(res(i,j+jone,k)+res(i,j-jone,k)) &
360 + b(3)*(res(i,j,k+kone)+res(i,j,k-kone)) &
365 call zero_boundaries (out)
370 end subroutine precond
376 SUBROUTINE sor (self, patch, phi1, d, d0, detailed_timer)
379 real(kind=KindScalarVar),
dimension(:,:,:),
pointer:: phi1
380 real(kind=KindScalarVar),
dimension(:,:,:),
pointer,
intent(in):: d
382 logical :: detailed_timer
384 integer :: i, j, k, sweep, ss, ione, jone, kone, m(3), iter, niter
385 real :: error, error0
386 real(8) :: a, b, b1, b2, b3, res, numer, denom, wc, usppt
387 integer,
save :: itimer=0
388 real(8),
dimension(:,:,:),
allocatable:: phi, source
391 if (detailed_timer)
call trace%begin (
'poisson_t%sor_iter', itimer=itimer)
393 associate(n=>patch%n,ds=>patch%ds,l=>patch%li,u=>patch%ui,gn=>patch%gn)
394 allocate (phi(gn(1),gn(2),gn(3)), source(gn(1),gn(2),gn(3)))
411 if (patch%n(i) > 1)
then 412 numer = numer + cos(pi / patch%n(i))
416 self%rjac2 = (numer / denom)**2
422 self%omega_sor = 2d0 / (1d0 + sqrt(1d0 - self%rjac2))
425 if (patch%is_periodic())
call patch%make_periodic (phi)
429 if (n(1) <= 1) ione = 0
430 if (n(2) <= 1) jone = 0
431 if (n(3) <= 1) kone = 0
435 source = fourpig * (d - d0)
438 if (n(1) > 1) a = a + 2. / ds(1)**2
439 if (n(2) > 1) a = a + 2. / ds(2)**2
440 if (n(3) > 1) a = a + 2. / ds(3)**2
442 b1 = a / ds(1)**2 * ione
443 b2 = a / ds(2)**2 * jone
444 b3 = a / ds(3)**2 * kone
451 ss = modulo(j + k + sweep,2)
457 res = b1 * (phi(i+ione,j,k) + phi(i-ione,j,k)) &
458 + b2 * (phi(i,j+jone,k) + phi(i,j-jone,k)) &
459 + b3 * (phi(i,j,k+kone) + phi(i,j,k-kone)) &
460 - phi(i,j,k)- a * source(i,j,k)
461 phi(i,j,k) = phi(i,j,k) + self%omega_sor * res
462 error = max(abs(res / (d(i,j,k) + floor)), error)
468 self%omega_sor = 1.0 / (1.0 - self%rjac2 * 0.5)
470 self%omega_sor = 1.0 / (1.0 - self%rjac2 * 0.25 * self%omega_sor)
479 if (iter==1) error0 = error
480 if (verbose > 1) print *,
'poisson_t%sor: iter, error =', iter, error
481 if (error < tolerance)
exit 485 if (iter==max_iter+1)
then 486 print *,
'poisson_t%sor: no convergence, niter, error =', max_iter, error
487 else if (verbose > 0 .and. io%master)
then 488 usppt = 1e6 * (wallclock() - wc) / product(shape(phi))
489 print
'(a,i7,i5,3x,2e12.3,f8.2)', &
490 'poisson_t%sor: id, niter, initial/final error, us/pt =', &
491 patch%id, iter, error0, error, usppt
498 if (detailed_timer)
call trace%end (itimer)
504 SUBROUTINE laplace (ds, phi, source, res)
506 real,
dimension(:,:,:) :: phi, source, res
508 integer :: ix, iy, iz, m(3), ione, jone, kone
514 if (
size(phi,1) <= 1) ione = 0
515 if (
size(phi,2) <= 1) jone = 0
516 if (
size(phi,3) <= 1) kone = 0
526 res(ix,iy,iz) = c(1)*(phi(ix+ione,iy,iz)+phi(ix-ione,iy,iz)-2d0*phi(ix,iy,iz)) &
527 + c(2)*(phi(ix,iy+jone,iz)+phi(ix,iy-jone,iz)-2d0*phi(ix,iy,iz)) &
528 + c(3)*(phi(ix,iy,iz+kone)+phi(ix,iy,iz-kone)-2d0*phi(ix,iy,iz)) &
533 END SUBROUTINE laplace
Support tic/toc timing, as in MATLAB, and accurate wallclock() function. The timing is generally much...
Define code units, in terms of (here) CGS units.
Template module for patches, which adds pointers to memory and mesh, and number of dimensions and var...
Methods for solving the Poisson equation; partially derived from routines developed by Troels Haugbøll...
Template module for mesh.