|
| 1 | +pure module function laplacian(rhs) result(laplacian_rhs) |
| 2 | + implicit none |
| 3 | + class(subdomain_t), intent(in) :: rhs |
| 4 | + type(subdomain_t) :: laplacian_rhs |
| 5 | + |
| 6 | + ! Local variables for array bounds and grid spacing |
| 7 | + integer :: my_nx, ny, nz, me, num_subdomains |
| 8 | + real :: dx_val, dy_val, dz_val, dx2_inv, dy2_inv, dz2_inv |
| 9 | + |
| 10 | + ! Get array dimensions and parallel image information |
| 11 | + my_nx = size(rhs%s_, 1) |
| 12 | + ny = size(rhs%s_, 2) |
| 13 | + nz = size(rhs%s_, 3) |
| 14 | + me = this_image() |
| 15 | + num_subdomains = num_images() |
| 16 | + |
| 17 | + ! Get grid spacing values |
| 18 | + dx_val = rhs%dx() |
| 19 | + dy_val = rhs%dy() |
| 20 | + dz_val = rhs%dz() |
| 21 | + |
| 22 | + ! Compute inverse squared grid spacings for efficiency |
| 23 | + dx2_inv = 1.0 / (dx_val * dx_val) |
| 24 | + dy2_inv = 1.0 / (dy_val * dy_val) |
| 25 | + dz2_inv = 1.0 / (dz_val * dz_val) |
| 26 | + |
| 27 | + ! Allocate result array with same shape as input |
| 28 | + allocate(laplacian_rhs%s_(my_nx, ny, nz)) |
| 29 | + |
| 30 | + ! Initialize result to zero |
| 31 | + laplacian_rhs%s_ = 0.0 |
| 32 | + |
| 33 | + ! Compute 2nd-order central difference Laplacian in interior points |
| 34 | + ! Laplacian = d²/dx² + d²/dy² + d²/dz² |
| 35 | + do concurrent (i = 2:my_nx-1, j = 2:ny-1, k = 2:nz-1) default(none) & |
| 36 | + shared(laplacian_rhs, rhs, dx2_inv, dy2_inv, dz2_inv, my_nx, ny, nz, me, num_subdomains) |
| 37 | + |
| 38 | + ! Second derivative in x-direction: (f(i+1) - 2*f(i) + f(i-1))/dx² |
| 39 | + ! Second derivative in y-direction: (f(j+1) - 2*f(j) + f(j-1))/dy² |
| 40 | + ! Second derivative in z-direction: (f(k+1) - 2*f(k) + f(k-1))/dz² |
| 41 | + laplacian_rhs%s_(i,j,k) = & |
| 42 | + dx2_inv * (rhs%s_(i+1,j,k) - 2.0*rhs%s_(i,j,k) + rhs%s_(i-1,j,k)) + & |
| 43 | + dy2_inv * (rhs%s_(i,j+1,k) - 2.0*rhs%s_(i,j,k) + rhs%s_(i,j-1,k)) + & |
| 44 | + dz2_inv * (rhs%s_(i,j,k+1) - 2.0*rhs%s_(i,j,k) + rhs%s_(i,j,k-1)) |
| 45 | + end do |
| 46 | + |
| 47 | + ! Handle boundary conditions at subdomain interfaces for x-direction |
| 48 | + ! Left boundary (i=1) - only compute if not the first subdomain |
| 49 | + if (me > 1) then |
| 50 | + do concurrent (j = 2:ny-1, k = 2:nz-1) default(none) & |
| 51 | + shared(laplacian_rhs, rhs, dx2_inv, dy2_inv, dz2_inv, ny, nz, halo_x) |
| 52 | + |
| 53 | + laplacian_rhs%s_(1,j,k) = & |
| 54 | + dx2_inv * (rhs%s_(2,j,k) - 2.0*rhs%s_(1,j,k) + halo_x(west,j,k)) + & |
| 55 | + dy2_inv * (rhs%s_(1,j+1,k) - 2.0*rhs%s_(1,j,k) + rhs%s_(1,j-1,k)) + & |
| 56 | + dz2_inv * (rhs%s_(1,j,k+1) - 2.0*rhs%s_(1,j,k) + rhs%s_(1,j,k-1)) |
| 57 | + end do |
| 58 | + end if |
| 59 | + |
| 60 | + ! Right boundary (i=my_nx) - only compute if not the last subdomain |
| 61 | + if (me < num_subdomains) then |
| 62 | + do concurrent (j = 2:ny-1, k = 2:nz-1) default(none) & |
| 63 | + shared(laplacian_rhs, rhs, dx2_inv, dy2_inv, dz2_inv, my_nx, ny, nz, halo_x, num_subdomains) |
| 64 | + |
| 65 | + laplacian_rhs%s_(my_nx,j,k) = & |
| 66 | + dx2_inv * (halo_x(east,j,k) - 2.0*rhs%s_(my_nx,j,k) + rhs%s_(my_nx-1,j,k)) + & |
| 67 | + dy2_inv * (rhs%s_(my_nx,j+1,k) - 2.0*rhs%s_(my_nx,j,k) + rhs%s_(my_nx,j-1,k)) + & |
| 68 | + dz2_inv * (rhs%s_(my_nx,j,k+1) - 2.0*rhs%s_(my_nx,j,k) + rhs%s_(my_nx,j,k-1)) |
| 69 | + end do |
| 70 | + end if |
| 71 | + |
| 72 | + ! Set boundary values to zero as specified |
| 73 | + ! Y boundaries |
| 74 | + laplacian_rhs%s_(:, 1, :) = 0.0 |
| 75 | + laplacian_rhs%s_(:, ny, :) = 0.0 |
| 76 | + |
| 77 | + ! Z boundaries |
| 78 | + laplacian_rhs%s_(:, :, 1) = 0.0 |
| 79 | + laplacian_rhs%s_(:, :, nz) = 0.0 |
| 80 | + |
| 81 | + ! X boundaries for first and last subdomains |
| 82 | + if (me == 1) then |
| 83 | + laplacian_rhs%s_(1, :, :) = 0.0 |
| 84 | + end if |
| 85 | + |
| 86 | + if (me == num_subdomains) then |
| 87 | + laplacian_rhs%s_(my_nx, :, :) = 0.0 |
| 88 | + end if |
| 89 | + |
| 90 | +end function laplacian |
0 commit comments