From 4d7f971466ac8c89d4cb43d4fb455485747e6d92 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Rui=20Ap=C3=B3stolo?= Date: Tue, 24 Oct 2023 17:41:20 +0100 Subject: [PATCH 01/62] temp --- src/3d/parcels/parcel_ellipsoid.f90 | 40 +++++++++++++++-------------- src/3d/parcels/parcel_merge.f90 | 13 +++++----- src/3d/stepper/rk_utils.f90 | 14 +++++----- 3 files changed, 36 insertions(+), 31 deletions(-) diff --git a/src/3d/parcels/parcel_ellipsoid.f90 b/src/3d/parcels/parcel_ellipsoid.f90 index 03ec2f208..573ce1836 100644 --- a/src/3d/parcels/parcel_ellipsoid.f90 +++ b/src/3d/parcels/parcel_ellipsoid.f90 @@ -35,7 +35,8 @@ module parcel_ellipsoid , I_B12 = 2 & ! index for B12 matrix component , I_B13 = 3 & ! index for B13 matrix component , I_B22 = 4 & ! index for B22 matrix component - , I_B23 = 5 ! index for B23 matrix component + , I_B23 = 5 & ! index for B23 matrix component + , I_B33 = 6 ! index for B33 matrix component integer :: IDX_ELL_VETA, IDX_ELL_VTAU @@ -132,9 +133,10 @@ end subroutine parcel_ellipsoid_deallocate ! @param[in] B = (B11, B12, B13, B22, B23) ! @param[in] volume of the parcel ! @returns the upper trinagular matrix - function get_full_matrix(B, volume) result(U) - double precision, intent(in) :: B(5) - double precision, intent(in) :: volume + function get_full_matrix(B) result(U) + ! function get_full_matrix(B, volume) result(U) + double precision, intent(in) :: B(6) + ! double precision, intent(in) :: volume double precision :: U(n_dim, n_dim) U(1, 1) = B(I_B11) @@ -145,7 +147,7 @@ function get_full_matrix(B, volume) result(U) U(2, 3) = B(I_B23) U(3, 1) = U(1, 3) U(3, 2) = U(2, 3) - U(3, 3) = get_B33(B, volume) + U(3, 3) = B(I_B33) end function get_full_matrix !:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: @@ -154,13 +156,13 @@ end function get_full_matrix ! @param[in] B = (B11, B12, B13, B22, B23) ! @param[in] volume of the parcel ! @returns all eigenvalues (sorted in descending order) - function get_eigenvalues(B, volume) result(D) + function get_eigenvalues(B) result(D) double precision, intent(in) :: B(5) - double precision, intent(in) :: volume + ! double precision, intent(in) :: volume double precision :: U(n_dim, n_dim) double precision :: D(n_dim) - U = get_full_matrix(B, volume) + U = get_full_matrix(B) call scherzinger_eigenvalues(U, D) @@ -178,12 +180,12 @@ end function get_eigenvalues function get_determinant(B, volume) result(detB) double precision, intent(in) :: B(5) double precision, intent(in) :: volume - double precision :: detB, B33 + double precision :: detB - B33 = get_B33(B, volume) + ! B33 = get_B33(B, volume) - detB = B(I_B11) * (B(I_B22) * B33 - B(I_B23) ** 2) & - - B(I_B12) * (B(I_B12) * B33 - B(I_B13) * B(I_B23)) & + detB = B(I_B11) * (B(I_B22) * B(I_B33) - B(I_B23) ** 2) & + - B(I_B12) * (B(I_B12) * B(I_B33) - B(I_B13) * B(I_B23)) & + B(I_B13) * (B(I_B12) * B(I_B23) - B(I_B13) * B(I_B22)) end function get_determinant @@ -195,12 +197,12 @@ end function get_determinant ! @param[in] B = (B11, B12, B13, B22, B23) ! @param[in] volume of the parcel ! @returns the eigenvectors - function get_eigenvectors(B, volume) result(V) - double precision, intent(in) :: B(5) - double precision, intent(in) :: volume + function get_eigenvectors(B) result(V) + double precision, intent(in) :: B(6) + ! double precision, intent(in) :: volume double precision :: U(n_dim, n_dim), D(n_dim), V(n_dim, n_dim) - U = get_full_matrix(B, volume) + U = get_full_matrix(B) call scherzinger_diagonalise(U, D, V) @@ -224,13 +226,13 @@ end function get_eigenvectors ! @param[in] volume of the parcel ! @param[out] D eigenvalues (sorted in descending order) ! @param[out] V eigenvectors - subroutine diagonalise(B, volume, D, V) - double precision, intent(in) :: B(5) + subroutine diagonalise(B, D, V) + double precision, intent(in) :: B(6) double precision, intent(in) :: volume double precision, intent(out) :: D(n_dim), V(n_dim, n_dim) double precision :: U(n_dim, n_dim) - U = get_full_matrix(B, volume) + U = get_full_matrix(B) call scherzinger_diagonalise(U, D, V) diff --git a/src/3d/parcels/parcel_merge.f90 b/src/3d/parcels/parcel_merge.f90 index 439537712..3820130e6 100644 --- a/src/3d/parcels/parcel_merge.f90 +++ b/src/3d/parcels/parcel_merge.f90 @@ -12,7 +12,8 @@ module parcel_merging , get_delx_across_periodic & , get_dely_across_periodic & , parcel_delete - use parcel_ellipsoid, only : get_B33, get_abc + use parcel_ellipsoid, only : get_abc + ! use parcel_ellipsoid, only : get_B33, get_abc use options, only : parcel #if defined (ENABLE_VERBOSE) && !defined (NDEBUG) use options, only : verbose @@ -104,7 +105,7 @@ subroutine do_group_merge(isma, iclo, n_merge, Bm, vm) integer :: loca(n_parcels) double precision :: x0(n_merge), y0(n_merge) double precision :: posm(3, n_merge) - double precision :: delx, vmerge, dely, delz, B33, mu + double precision :: delx, vmerge, dely, delz, mu double precision :: buoym(n_merge), vortm(3, n_merge) #ifndef ENABLE_DRY_MODE double precision :: hum(n_merge) @@ -222,7 +223,7 @@ subroutine do_group_merge(isma, iclo, n_merge, Bm, vm) vmerge = one / vm(l) - B33 = get_B33(parcels%B(:, ic), parcels%volume(ic)) + ! B33 = get_B33(parcels%B(:, ic), parcels%volume(ic)) delx = get_delx_across_periodic(parcels%position(1, ic), posm(1, l)) dely = get_dely_across_periodic(parcels%position(2, ic), posm(2, l)) @@ -236,7 +237,7 @@ subroutine do_group_merge(isma, iclo, n_merge, Bm, vm) Bm(3, l) = mu * (five * delx * delz + parcels%B(3, ic)) Bm(4, l) = mu * (five * dely ** 2 + parcels%B(4, ic)) Bm(5, l) = mu * (five * dely * delz + parcels%B(5, ic)) - Bm(6, l) = mu * (five * delz ** 2 + B33) + Bm(6, l) = mu * (five * delz ** 2 + parcels%B(6, ic)) parcels%volume(ic) = vm(l) parcels%position(1, ic) = posm(1, l) @@ -260,7 +261,7 @@ subroutine do_group_merge(isma, iclo, n_merge, Bm, vm) dely = get_dely_across_periodic(parcels%position(2, is), posm(2, n)) delz = parcels%position(3, is) - posm(3, n) - B33 = get_B33(parcels%B(:, is), parcels%volume(is)) + ! B33 = get_B33(parcels%B(:, is), parcels%volume(is)) ! volume fraction V_{is} / V mu = vmerge * parcels%volume(is) @@ -270,7 +271,7 @@ subroutine do_group_merge(isma, iclo, n_merge, Bm, vm) Bm(3, n) = Bm(3, n) + mu * (five * delx * delz + parcels%B(3, is)) Bm(4, n) = Bm(4, n) + mu * (five * dely ** 2 + parcels%B(4, is)) Bm(5, n) = Bm(5, n) + mu * (five * dely * delz + parcels%B(5, is)) - Bm(6, n) = Bm(6, n) + mu * (five * delz ** 2 + B33) + Bm(6, n) = Bm(6, n) + mu * (five * delz ** 2 + parcels%B(5, is)) enddo end subroutine do_group_merge diff --git a/src/3d/stepper/rk_utils.f90 b/src/3d/stepper/rk_utils.f90 index 17b50d4da..0fb6551f1 100644 --- a/src/3d/stepper/rk_utils.f90 +++ b/src/3d/stepper/rk_utils.f90 @@ -1,6 +1,6 @@ module rk_utils use dimensions, only : n_dim, I_X, I_Y, I_Z - use parcel_ellipsoid, only : get_B33, I_B11, I_B12, I_B13, I_B22, I_B23 + use parcel_ellipsoid, only : get_B33, I_B11, I_B12, I_B13, I_B22, I_B23, I_B33 use fields, only : velgradg, tbuoyg, vortg, I_DUDX, I_DUDY, I_DVDY, I_DWDX, I_DWDY use field_mpi, only : field_halo_fill use constants, only : zero, one, two, f12 @@ -27,11 +27,11 @@ module rk_utils ! @param[in] volume is the parcel volume ! @returns dB/dt in Bout function get_dBdt(Bin, S, vorticity, volume) result(Bout) - double precision, intent(in) :: Bin(I_B23) + double precision, intent(in) :: Bin(I_B33) double precision, intent(in) :: S(5) double precision, intent(in) :: vorticity(n_dim) double precision, intent(in) :: volume - double precision :: Bout(5), B33 + double precision :: Bout(I_B33) double precision :: dudz, dvdx, dvdz, dwdz ! du/dz = \eta + dw/dx @@ -46,7 +46,7 @@ function get_dBdt(Bin, S, vorticity, volume) result(Bout) ! dw/dz = - (du/dx + dv/dy) dwdz = - (S(I_DUDX) + S(I_DVDY)) - B33 = get_B33(Bin, volume) + ! B33 = get_B33(Bin, volume) ! dB11/dt = 2 * (du/dx * B11 + du/dy * B12 + du/dz * B13) Bout(I_B11) = two * (S(I_DUDX) * Bin(I_B11) + S(I_DUDY) * Bin(I_B12) + dudz * Bin(I_B13)) @@ -63,7 +63,7 @@ function get_dBdt(Bin, S, vorticity, volume) result(Bout) + S(I_DWDY) * Bin(I_B12) & ! + dw/dy * B12 - S(I_DVDY) * Bin(I_B13) & ! - dv/dy * B13 + S(I_DUDY) * Bin(I_B23) & ! + du/dy * B23 - + dudz * B33 ! + du/dz * B33 + + dudz * Bin(I_B33) ! + du/dz * B33 ! dB22/dt = 2 * (dv/dx * B12 + dv/dy * B22 + dv/dz * B23) Bout(I_B22) = two * (dvdx * Bin(I_B12) + S(I_DVDY) * Bin(I_B22) + dvdz * Bin(I_B23)) @@ -73,7 +73,9 @@ function get_dBdt(Bin, S, vorticity, volume) result(Bout) + dvdx * Bin(I_B13) & ! + dv/dx * B13 + S(I_DWDY) * Bin(I_B22) & ! + dw/dy * B22 - S(I_DUDX) * Bin(I_B23) & ! - du/dx * B23 - + dvdz * B33 ! + dv/dz * B33 + + dvdz * Bin(I_B33) ! + dv/dz * B33 + + Bout(I_B33) = get_B33(Bout(I_B11:I_B23), volume) end function get_dBdt ! Estimate a suitable time step based on the velocity strain From e2e39ec874d526906a45df3cd7c7a928655c0991 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Rui=20Ap=C3=B3stolo?= Date: Mon, 30 Oct 2023 08:54:55 +0000 Subject: [PATCH 02/62] Changed size of B matrix variables --- src/3d/parcels/parcel_ellipsoid.f90 | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/src/3d/parcels/parcel_ellipsoid.f90 b/src/3d/parcels/parcel_ellipsoid.f90 index 573ce1836..4c3d47907 100644 --- a/src/3d/parcels/parcel_ellipsoid.f90 +++ b/src/3d/parcels/parcel_ellipsoid.f90 @@ -153,11 +153,11 @@ end function get_full_matrix !:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: ! Obtain all eigenvalues sorted in descending order - ! @param[in] B = (B11, B12, B13, B22, B23) + ! @param[in] B = (B11, B12, B13, B22, B23, B33) ! @param[in] volume of the parcel ! @returns all eigenvalues (sorted in descending order) function get_eigenvalues(B) result(D) - double precision, intent(in) :: B(5) + double precision, intent(in) :: B(6) ! double precision, intent(in) :: volume double precision :: U(n_dim, n_dim) double precision :: D(n_dim) @@ -178,7 +178,7 @@ end function get_eigenvalues !:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: function get_determinant(B, volume) result(detB) - double precision, intent(in) :: B(5) + double precision, intent(in) :: B(6) double precision, intent(in) :: volume double precision :: detB @@ -306,7 +306,7 @@ end function get_aspect_ratio function get_ellipsoid_points(position, volume, B, n, l_reuse) result(points) double precision, intent(in) :: position(n_dim) double precision, intent(in) :: volume - double precision, intent(in) :: B(5) ! B11, B12, B13, B22, B23 + double precision, intent(in) :: B(6) ! B11, B12, B13, B22, B23 integer, optional, intent(in) :: n logical, optional, intent(in) :: l_reuse double precision :: Veta(n_dim), Vtau(n_dim), D(n_dim), V(n_dim, n_dim) @@ -354,7 +354,7 @@ end function get_ellipsoid_points !:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: function get_angles(B, volume) result(angles) - double precision, intent(in) :: B(5) + double precision, intent(in) :: B(6) double precision, intent(in) :: volume double precision :: evec(n_dim, n_dim) double precision :: angles(2) ! (/azimuth, polar/) From 913ace008352b11dae47b66c1de918126aa04714 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Rui=20Ap=C3=B3stolo?= Date: Mon, 30 Oct 2023 09:01:59 +0000 Subject: [PATCH 03/62] Removed unused volume parameters in function calls --- src/3d/parcels/parcel_ellipsoid.f90 | 20 ++++++++++---------- src/3d/parcels/parcel_interpl.f90 | 6 +++--- 2 files changed, 13 insertions(+), 13 deletions(-) diff --git a/src/3d/parcels/parcel_ellipsoid.f90 b/src/3d/parcels/parcel_ellipsoid.f90 index 4c3d47907..6096ef99b 100644 --- a/src/3d/parcels/parcel_ellipsoid.f90 +++ b/src/3d/parcels/parcel_ellipsoid.f90 @@ -177,9 +177,9 @@ end function get_eigenvalues !:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: - function get_determinant(B, volume) result(detB) + function get_determinant(B) result(detB) double precision, intent(in) :: B(6) - double precision, intent(in) :: volume + ! double precision, intent(in) :: volume double precision :: detB ! B33 = get_B33(B, volume) @@ -228,7 +228,7 @@ end function get_eigenvectors ! @param[out] V eigenvectors subroutine diagonalise(B, D, V) double precision, intent(in) :: B(6) - double precision, intent(in) :: volume + ! double precision, intent(in) :: volume double precision, intent(out) :: D(n_dim), V(n_dim, n_dim) double precision :: U(n_dim, n_dim) @@ -303,9 +303,9 @@ end function get_aspect_ratio ! @param[in] volume of the parcel ! @param[in] B matrix elements of the parcel ! @returns the parcel support points - function get_ellipsoid_points(position, volume, B, n, l_reuse) result(points) + function get_ellipsoid_points(position, B, n, l_reuse) result(points) double precision, intent(in) :: position(n_dim) - double precision, intent(in) :: volume + ! double precision, intent(in) :: volume double precision, intent(in) :: B(6) ! B11, B12, B13, B22, B23 integer, optional, intent(in) :: n logical, optional, intent(in) :: l_reuse @@ -319,7 +319,7 @@ function get_ellipsoid_points(position, volume, B, n, l_reuse) result(points) Veta = Vetas(:, n) Vtau = Vtaus(:, n) else - call diagonalise(B, volume, D, V) + call diagonalise(B, D, V) Veta = dsqrt(dabs(D(I_X) - D(I_Z))) * rho * V(:, I_X) Vtau = dsqrt(dabs(D(I_Y) - D(I_Z))) * rho * V(:, I_Y) @@ -329,7 +329,7 @@ function get_ellipsoid_points(position, volume, B, n, l_reuse) result(points) else ! (/a2, b2, c2/) with a >= b >= c ! D = (/a2, b2, c2/) - call diagonalise(B, volume, D, V) + call diagonalise(B, D, V) Veta = dsqrt(dabs(D(I_X) - D(I_Z))) * rho * V(:, I_X) Vtau = dsqrt(dabs(D(I_Y) - D(I_Z))) * rho * V(:, I_Y) @@ -353,13 +353,13 @@ end function get_ellipsoid_points !:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: - function get_angles(B, volume) result(angles) + function get_angles(B) result(angles) double precision, intent(in) :: B(6) - double precision, intent(in) :: volume + ! double precision, intent(in) :: volume double precision :: evec(n_dim, n_dim) double precision :: angles(2) ! (/azimuth, polar/) - evec = get_eigenvectors(B, volume) + evec = get_eigenvectors(B) ! azimuthal angle angles(I_X) = datan2(evec(I_Y, I_X), evec(I_X, I_X)) diff --git a/src/3d/parcels/parcel_interpl.f90 b/src/3d/parcels/parcel_interpl.f90 index 67529c80b..0b9082297 100644 --- a/src/3d/parcels/parcel_interpl.f90 +++ b/src/3d/parcels/parcel_interpl.f90 @@ -88,7 +88,7 @@ subroutine vol2grid(l_reuse) pvol = parcels%volume(n) points = get_ellipsoid_points(parcels%position(:, n), & - pvol, parcels%B(:, n), & + parcels%B(:, n), & n, l_reuse) ! we have 4 points per ellipsoid @@ -179,7 +179,7 @@ subroutine par2grid(l_reuse) btot = btot - bfsq * parcels%position(3, n) #endif points = get_ellipsoid_points(parcels%position(:, n), & - pvol, parcels%B(:, n), n, l_reuse) + parcels%B(:, n), n, l_reuse) call get_index(parcels%position(:, n), i, j, k) nparg(k, j, i) = nparg(k, j, i) + 1 @@ -422,7 +422,7 @@ subroutine grid2par(add) parcels%strain(:, n) = zero points = get_ellipsoid_points(parcels%position(:, n), & - parcels%volume(n), parcels%B(:, n), n) + parcels%B(:, n), n) do p = 1, 4 ! get interpolation weights and mesh indices From cce6b591c3865ae487e9d94d65902e2f879d58b9 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Rui=20Ap=C3=B3stolo?= Date: Mon, 30 Oct 2023 09:06:57 +0000 Subject: [PATCH 04/62] Changes to parcel_split --- src/3d/parcels/parcel_split.f90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/3d/parcels/parcel_split.f90 b/src/3d/parcels/parcel_split.f90 index 6be004582..835e2d8d7 100644 --- a/src/3d/parcels/parcel_split.f90 +++ b/src/3d/parcels/parcel_split.f90 @@ -38,7 +38,7 @@ module parcel_split_mod ! Split elongated parcels (semi-major axis larger than amax) or ! parcels with aspect ratios larger than parcel%lambda_max. subroutine parcel_split - double precision :: B(5) + double precision :: B(6) double precision :: vol, lam double precision :: D(3), V(3, 3) integer :: last_index, n_indices @@ -117,7 +117,7 @@ subroutine parcel_split B = parcels%B(:, n) vol = parcels%volume(n) - call diagonalise(B, vol, D, V) + call diagonalise(B, D, V) pid(n) = 0 From 3b2796314542008809c494370d5aa38bf1f29d6a Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Rui=20Ap=C3=B3stolo?= Date: Mon, 30 Oct 2023 09:08:30 +0000 Subject: [PATCH 05/62] Changes to parcel_split --- src/3d/parcels/parcel_split.f90 | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/3d/parcels/parcel_split.f90 b/src/3d/parcels/parcel_split.f90 index 835e2d8d7..bf4d03657 100644 --- a/src/3d/parcels/parcel_split.f90 +++ b/src/3d/parcels/parcel_split.f90 @@ -39,7 +39,7 @@ module parcel_split_mod ! parcels with aspect ratios larger than parcel%lambda_max. subroutine parcel_split double precision :: B(6) - double precision :: vol, lam + double precision :: lam double precision :: D(3), V(3, 3) integer :: last_index, n_indices integer :: grown_size, shrunk_size, n_required @@ -56,12 +56,12 @@ subroutine parcel_split !------------------------------------------------------------------ ! Check which parcels split and store the indices in *pid*: !$omp parallel default(shared) - !$omp do private(n, B, vol, lam, D) + !$omp do private(n, B, lam, D) do n = 1, n_parcels B = parcels%B(:, n) - vol = parcels%volume(n) + ! vol = parcels%volume(n) - D = get_eigenvalues(B, vol) + D = get_eigenvalues(B) ! evaluate maximum aspect ratio (a2 >= b2 >= c2) lam = get_aspect_ratio(D) From ea2956f66562263185bd7de3263722695d353b01 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Rui=20Ap=C3=B3stolo?= Date: Mon, 30 Oct 2023 09:09:07 +0000 Subject: [PATCH 06/62] Changes to parcel_split --- src/3d/parcels/parcel_split.f90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/3d/parcels/parcel_split.f90 b/src/3d/parcels/parcel_split.f90 index bf4d03657..b2f9ba178 100644 --- a/src/3d/parcels/parcel_split.f90 +++ b/src/3d/parcels/parcel_split.f90 @@ -39,7 +39,7 @@ module parcel_split_mod ! parcels with aspect ratios larger than parcel%lambda_max. subroutine parcel_split double precision :: B(6) - double precision :: lam + double precision :: vol, lam double precision :: D(3), V(3, 3) integer :: last_index, n_indices integer :: grown_size, shrunk_size, n_required From f4df232de29cbbb5ef6a44621b9d41fc9dac508d Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Rui=20Ap=C3=B3stolo?= Date: Mon, 30 Oct 2023 09:10:15 +0000 Subject: [PATCH 07/62] Changes to parcel_diagnostics --- src/3d/parcels/parcel_diagnostics.f90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/3d/parcels/parcel_diagnostics.f90 b/src/3d/parcels/parcel_diagnostics.f90 index ee5aabb67..559d395dd 100644 --- a/src/3d/parcels/parcel_diagnostics.f90 +++ b/src/3d/parcels/parcel_diagnostics.f90 @@ -84,7 +84,7 @@ subroutine calculate_parcel_diagnostics parcel_stats(IDX_ENSTROPHY) = parcel_stats(IDX_ENSTROPHY) & + (vor(1) ** 2 + vor(2) ** 2 + vor(3) ** 2) * vol - evals = get_eigenvalues(parcels%B(:, n), parcels%volume(n)) + evals = get_eigenvalues(parcels%B(:, n)) lam = get_aspect_ratio(evals) parcel_stats(IDX_AVG_LAM) = parcel_stats(IDX_AVG_LAM) + lam @@ -103,7 +103,7 @@ subroutine calculate_parcel_diagnostics #ifndef NDEBUG !$omp critical - if (abs(get_determinant(parcels%B(:, n), vol) / get_abc(vol) ** 2 - one) > thres) then + if (abs(get_determinant(parcels%B(:, n)) / get_abc(vol) ** 2 - one) > thres) then call mpi_exit_on_error("Parcel determinant not preserved!") endif !$omp end critical From 6f0bee5bac9b00d2f1204a48ca97e1ea4ad013be Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Rui=20Ap=C3=B3stolo?= Date: Mon, 30 Oct 2023 09:11:28 +0000 Subject: [PATCH 08/62] Changes to parcel_init --- src/3d/parcels/parcel_init.f90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/3d/parcels/parcel_init.f90 b/src/3d/parcels/parcel_init.f90 index 837e95cf3..ae66cf325 100644 --- a/src/3d/parcels/parcel_init.f90 +++ b/src/3d/parcels/parcel_init.f90 @@ -163,7 +163,7 @@ subroutine init_refine(lam) ! do refining by splitting do while (lam >= parcel%lambda_max) call parcel_split - evals = get_eigenvalues(parcels%B(1, :), parcels%volume(1)) + evals = get_eigenvalues(parcels%B(1, :)) lam = dsqrt(evals(1) / evals(3)) end do end subroutine init_refine From f7f59928e99a0b054ca8bffb63aabd1018bd675f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Rui=20Ap=C3=B3stolo?= Date: Mon, 30 Oct 2023 09:30:15 +0000 Subject: [PATCH 09/62] Change to parcel container --- src/3d/parcels/parcel_container.f90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/3d/parcels/parcel_container.f90 b/src/3d/parcels/parcel_container.f90 index faa92f3e9..e4f1fc4f5 100644 --- a/src/3d/parcels/parcel_container.f90 +++ b/src/3d/parcels/parcel_container.f90 @@ -304,7 +304,7 @@ subroutine parcel_alloc(num) allocate(parcels%position(3, num)) allocate(parcels%vorticity(3, num)) - allocate(parcels%B(5, num)) + allocate(parcels%B(6, num)) allocate(parcels%volume(num)) allocate(parcels%buoyancy(num)) #ifndef ENABLE_DRY_MODE From 4dceb95dfee0c4fd0a5b82fb14d1165c4b09588d Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Rui=20Ap=C3=B3stolo?= Date: Mon, 30 Oct 2023 09:30:33 +0000 Subject: [PATCH 10/62] Change to test_mpi_parcel_split --- unit-tests/3d/test_mpi_parcel_split.f90 | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/unit-tests/3d/test_mpi_parcel_split.f90 b/unit-tests/3d/test_mpi_parcel_split.f90 index e8766ec39..0fb0914ba 100644 --- a/unit-tests/3d/test_mpi_parcel_split.f90 +++ b/unit-tests/3d/test_mpi_parcel_split.f90 @@ -141,7 +141,7 @@ program test_mpi_parcel_split subroutine fill_shape(m, angle) integer, intent(in) :: m double precision, intent(in) :: angle - double precision :: B11, B12, B13, B22, B23, st, ct, sp, cp + double precision :: B11, B12, B13, B22, B23, B33, st, ct, sp, cp st = dsin(angle) ct = dcos(angle) @@ -159,6 +159,7 @@ subroutine fill_shape(m, angle) parcels%B(3, m) = B13 parcels%B(4, m) = B22 parcels%B(5, m) = B23 + parcels%B(6, m) = B33 end subroutine fill_shape end program test_mpi_parcel_split From 0993089353b600f54fd8f110b73f430d72564dc3 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Rui=20Ap=C3=B3stolo?= Date: Mon, 30 Oct 2023 09:46:18 +0000 Subject: [PATCH 11/62] Change to parcel_bc --- src/3d/parcels/parcel_bc.f90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/3d/parcels/parcel_bc.f90 b/src/3d/parcels/parcel_bc.f90 index 770701ff0..45ba38f15 100644 --- a/src/3d/parcels/parcel_bc.f90 +++ b/src/3d/parcels/parcel_bc.f90 @@ -24,7 +24,7 @@ end subroutine apply_periodic_bc ! @param[inout] position vector of parcel ! @param[inout] B matrix of parcel pure subroutine apply_reflective_bc(position, B) - double precision, intent(inout) :: position(3), B(5) + double precision, intent(inout) :: position(3), B(6) if (position(3) > upper(3)) then position(3) = two * upper(3) - position(3) From 1b79fc5f6744172a7570af430d5ebe89cbc24f0e Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Rui=20Ap=C3=B3stolo?= Date: Mon, 30 Oct 2023 09:49:51 +0000 Subject: [PATCH 12/62] Change to test_mpi_parcel_split --- unit-tests/3d/test_mpi_parcel_split.f90 | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/unit-tests/3d/test_mpi_parcel_split.f90 b/unit-tests/3d/test_mpi_parcel_split.f90 index 0fb0914ba..e8766ec39 100644 --- a/unit-tests/3d/test_mpi_parcel_split.f90 +++ b/unit-tests/3d/test_mpi_parcel_split.f90 @@ -141,7 +141,7 @@ program test_mpi_parcel_split subroutine fill_shape(m, angle) integer, intent(in) :: m double precision, intent(in) :: angle - double precision :: B11, B12, B13, B22, B23, B33, st, ct, sp, cp + double precision :: B11, B12, B13, B22, B23, st, ct, sp, cp st = dsin(angle) ct = dcos(angle) @@ -159,7 +159,6 @@ subroutine fill_shape(m, angle) parcels%B(3, m) = B13 parcels%B(4, m) = B22 parcels%B(5, m) = B23 - parcels%B(6, m) = B33 end subroutine fill_shape end program test_mpi_parcel_split From fe71797dfb9d670b36a7a4e594e454e4faacf057 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Rui=20Ap=C3=B3stolo?= Date: Tue, 7 Nov 2023 11:45:39 +0000 Subject: [PATCH 13/62] Added B33 indexes to parcel communication in parcel_container --- src/3d/parcels/parcel_container.f90 | 28 ++++++++++++++++------------ 1 file changed, 16 insertions(+), 12 deletions(-) diff --git a/src/3d/parcels/parcel_container.f90 b/src/3d/parcels/parcel_container.f90 index e4f1fc4f5..f8ca33100 100644 --- a/src/3d/parcels/parcel_container.f90 +++ b/src/3d/parcels/parcel_container.f90 @@ -36,6 +36,7 @@ module parcel_container IDX_B13, & ! B13 shape matrix element IDX_B22, & ! B22 shape matrix element IDX_B23, & ! B23 shape matrix element + IDX_B33, & ! B33 shape matrix element IDX_VOL, & ! volume IDX_BUO, & ! buoyancy #ifndef ENABLE_DRY_MODE @@ -52,6 +53,7 @@ module parcel_container IDX_RK4_DB13, & ! RK4 variable for B13 IDX_RK4_DB22, & ! RK4 variable for B22 IDX_RK4_DB23, & ! RK4 variable for B23 + IDX_RK4_DB33, & ! RK4 variable for B33 IDX_RK4_DUDX, & ! RK4 variable du/dx IDX_RK4_DUDY, & ! RK4 variable du/dy IDX_RK4_DVDY, & ! RK4 variable dv/dy @@ -107,8 +109,9 @@ subroutine set_buffer_indices IDX_B13 = 9 ! B13 shape matrix element IDX_B22 = 10 ! B22 shape matrix element IDX_B23 = 11 ! B23 shape matrix element - IDX_VOL = 12 ! volume - IDX_BUO = 13 ! buoyancy + IDX_B33 = 12 ! B33 shape matrix element + IDX_VOL = 13 ! volume + IDX_BUO = 14 ! buoyancy i = IDX_BUO + 1 #ifndef ENABLE_DRY_MODE @@ -128,13 +131,14 @@ subroutine set_buffer_indices IDX_RK4_DB13 = i + 8 IDX_RK4_DB22 = i + 9 IDX_RK4_DB23 = i + 10 - IDX_RK4_DUDX = i + 11 - IDX_RK4_DUDY = i + 12 - IDX_RK4_DVDY = i + 13 - IDX_RK4_DWDX = i + 14 - IDX_RK4_DWDY = i + 15 + IDX_RK4_DB33 = i + 11 + IDX_RK4_DUDX = i + 12 + IDX_RK4_DUDY = i + 13 + IDX_RK4_DVDY = i + 14 + IDX_RK4_DWDX = i + 15 + IDX_RK4_DWDY = i + 16 - i = i + 16 + i = i + 17 n_par_attrib = set_ellipsoid_buffer_indices(i) @@ -359,7 +363,7 @@ subroutine parcel_serialize(n, buffer) buffer(IDX_X_POS:IDX_Z_POS) = parcels%position(:, n) buffer(IDX_X_VOR:IDX_Z_VOR) = parcels%vorticity(:, n) - buffer(IDX_B11:IDX_B23) = parcels%B(:, n) + buffer(IDX_B11:IDX_B33) = parcels%B(:, n) buffer(IDX_VOL) = parcels%volume(n) buffer(IDX_BUO) = parcels%buoyancy(n) #ifndef ENABLE_DRY_MODE @@ -368,7 +372,7 @@ subroutine parcel_serialize(n, buffer) ! LS-RK4 variables: buffer(IDX_RK4_X_DPOS:IDX_RK4_Z_DPOS) = parcels%delta_pos(:, n) buffer(IDX_RK4_X_DVOR:IDX_RK4_Z_DVOR) = parcels%delta_vor(:, n) - buffer(IDX_RK4_DB11:IDX_RK4_DB23) = parcels%delta_b(:, n) + buffer(IDX_RK4_DB11:IDX_RK4_DB33) = parcels%delta_b(:, n) buffer(IDX_RK4_DUDX:IDX_RK4_DWDY) = parcels%strain(:, n) call parcel_ellipsoid_serialize(n, buffer) @@ -383,7 +387,7 @@ subroutine parcel_deserialize(n, buffer) parcels%position(:, n) = buffer(IDX_X_POS:IDX_Z_POS) parcels%vorticity(:, n) = buffer(IDX_X_VOR:IDX_Z_VOR) - parcels%B(:, n) = buffer(IDX_B11:IDX_B23) + parcels%B(:, n) = buffer(IDX_B11:IDX_B33) parcels%volume(n) = buffer(IDX_VOL) parcels%buoyancy(n) = buffer(IDX_BUO) #ifndef ENABLE_DRY_MODE @@ -392,7 +396,7 @@ subroutine parcel_deserialize(n, buffer) ! LS-RK4 variables: parcels%delta_pos(:, n) = buffer(IDX_RK4_X_DPOS:IDX_RK4_Z_DPOS) parcels%delta_vor(:, n) = buffer(IDX_RK4_X_DVOR:IDX_RK4_Z_DVOR) - parcels%delta_b(:, n) = buffer(IDX_RK4_DB11:IDX_RK4_DB23) + parcels%delta_b(:, n) = buffer(IDX_RK4_DB11:IDX_RK4_DB33) parcels%strain(:, n) = buffer(IDX_RK4_DUDX:IDX_RK4_DWDY) call parcel_ellipsoid_deserialize(n, buffer) From a5036428d3bbb21054d1e4e7e60eb319747ada11 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Rui=20Ap=C3=B3stolo?= Date: Tue, 7 Nov 2023 12:18:28 +0000 Subject: [PATCH 14/62] Changed delta_b size --- src/3d/parcels/parcel_container.f90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/3d/parcels/parcel_container.f90 b/src/3d/parcels/parcel_container.f90 index f8ca33100..2302f10ab 100644 --- a/src/3d/parcels/parcel_container.f90 +++ b/src/3d/parcels/parcel_container.f90 @@ -70,7 +70,7 @@ module parcel_container vorticity, & B ! B matrix entries; ordering: ! B(:, 1) = B11, B(:, 2) = B12, B(:, 3) = B13 - ! B(:, 4) = B22, B(:, 5) = B23 + ! B(:, 4) = B22, B(:, 5) = B23, B(:, 6) = B33 double precision, allocatable, dimension(:) :: & volume, & @@ -320,7 +320,7 @@ subroutine parcel_alloc(num) allocate(parcels%delta_pos(3, num)) allocate(parcels%delta_vor(3, num)) allocate(parcels%strain(5, num)) - allocate(parcels%delta_b(5, num)) + allocate(parcels%delta_b(6, num)) end subroutine parcel_alloc From 5f7521e47358229da2a58842d5df8d2b94efdde8 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Rui=20Ap=C3=B3stolo?= Date: Tue, 7 Nov 2023 14:14:14 +0000 Subject: [PATCH 15/62] Docs --- src/3d/parcels/parcel_ellipsoid.f90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/3d/parcels/parcel_ellipsoid.f90 b/src/3d/parcels/parcel_ellipsoid.f90 index 6096ef99b..2991ea943 100644 --- a/src/3d/parcels/parcel_ellipsoid.f90 +++ b/src/3d/parcels/parcel_ellipsoid.f90 @@ -130,7 +130,7 @@ end subroutine parcel_ellipsoid_deallocate !:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: ! Obtain the parcel shape matrix. - ! @param[in] B = (B11, B12, B13, B22, B23) + ! @param[in] B = (B11, B12, B13, B22, B23, B33) ! @param[in] volume of the parcel ! @returns the upper trinagular matrix function get_full_matrix(B) result(U) @@ -306,7 +306,7 @@ end function get_aspect_ratio function get_ellipsoid_points(position, B, n, l_reuse) result(points) double precision, intent(in) :: position(n_dim) ! double precision, intent(in) :: volume - double precision, intent(in) :: B(6) ! B11, B12, B13, B22, B23 + double precision, intent(in) :: B(6) ! B11, B12, B13, B22, B23, B33 integer, optional, intent(in) :: n logical, optional, intent(in) :: l_reuse double precision :: Veta(n_dim), Vtau(n_dim), D(n_dim), V(n_dim, n_dim) From d7a559abd264ee645e87818753cde16c31d92148 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Rui=20Ap=C3=B3stolo?= Date: Wed, 8 Nov 2023 14:27:10 +0000 Subject: [PATCH 16/62] Modified unit test --- unit-tests/3d/test_mpi_gradient_correction_3d.f90 | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/unit-tests/3d/test_mpi_gradient_correction_3d.f90 b/unit-tests/3d/test_mpi_gradient_correction_3d.f90 index b8286e128..b4de335da 100644 --- a/unit-tests/3d/test_mpi_gradient_correction_3d.f90 +++ b/unit-tests/3d/test_mpi_gradient_correction_3d.f90 @@ -16,7 +16,7 @@ program test_mpi_gradient_correction_3d , vort_corr_timer & , init_parcel_correction use parcel_interpl, only : vol2grid, halo_swap_timer - use parcel_ellipsoid, only : get_abc + use parcel_ellipsoid, only : get_abc, get_b33 use parcel_init, only : init_regular_positions use parameters, only : lower, extent, update_parameters, vcell, nx, ny, nz, dx use fields, only : volg, field_default @@ -117,6 +117,9 @@ program test_mpi_gradient_correction_3d ! b22 parcels%B(4, :) = parcels%B(1, :) + ! b33 needs to be calculated here + parcels%B(6, :) = get_b33(parcels%B(0:5, :), parcels%volume) + call vol2grid volg(0:nz, lo(2):hi(2), lo(1):hi(1)) = abs(volg(0:nz, lo(2):hi(2), lo(1):hi(1)) / vcell - one) From ad987dff2a1c135d583cd422f276abd285672651 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Rui=20Ap=C3=B3stolo?= Date: Wed, 8 Nov 2023 14:32:14 +0000 Subject: [PATCH 17/62] Rank mismatch --- src/3d/parcels/parcel_ellipsoid.f90 | 2 +- unit-tests/3d/test_mpi_gradient_correction_3d.f90 | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/src/3d/parcels/parcel_ellipsoid.f90 b/src/3d/parcels/parcel_ellipsoid.f90 index 2991ea943..0c1d49b03 100644 --- a/src/3d/parcels/parcel_ellipsoid.f90 +++ b/src/3d/parcels/parcel_ellipsoid.f90 @@ -251,7 +251,7 @@ end subroutine diagonalise ! @param[in] B = (B11, B12, B13, B22, B23) ! @param[in] volume of the parcel ! @returns B33 - function get_B33(B, volume) result(B33) + pure elemental function get_B33(B, volume) result(B33) double precision, intent(in) :: B(5) double precision, intent(in) :: volume double precision :: abc diff --git a/unit-tests/3d/test_mpi_gradient_correction_3d.f90 b/unit-tests/3d/test_mpi_gradient_correction_3d.f90 index b4de335da..327422210 100644 --- a/unit-tests/3d/test_mpi_gradient_correction_3d.f90 +++ b/unit-tests/3d/test_mpi_gradient_correction_3d.f90 @@ -118,7 +118,7 @@ program test_mpi_gradient_correction_3d parcels%B(4, :) = parcels%B(1, :) ! b33 needs to be calculated here - parcels%B(6, :) = get_b33(parcels%B(0:5, :), parcels%volume) + parcels%B(6, :) = get_b33(parcels%B(0:5, :), parcels%volume(:)) call vol2grid From 9188c5423b0358acdb118d568f0f3daa982c3ac0 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Rui=20Ap=C3=B3stolo?= Date: Wed, 8 Nov 2023 14:49:02 +0000 Subject: [PATCH 18/62] get-B33 calls MPI, can't be pure function --- src/3d/parcels/parcel_ellipsoid.f90 | 2 +- unit-tests/3d/test_mpi_gradient_correction_3d.f90 | 4 +++- 2 files changed, 4 insertions(+), 2 deletions(-) diff --git a/src/3d/parcels/parcel_ellipsoid.f90 b/src/3d/parcels/parcel_ellipsoid.f90 index 0c1d49b03..2991ea943 100644 --- a/src/3d/parcels/parcel_ellipsoid.f90 +++ b/src/3d/parcels/parcel_ellipsoid.f90 @@ -251,7 +251,7 @@ end subroutine diagonalise ! @param[in] B = (B11, B12, B13, B22, B23) ! @param[in] volume of the parcel ! @returns B33 - pure elemental function get_B33(B, volume) result(B33) + function get_B33(B, volume) result(B33) double precision, intent(in) :: B(5) double precision, intent(in) :: volume double precision :: abc diff --git a/unit-tests/3d/test_mpi_gradient_correction_3d.f90 b/unit-tests/3d/test_mpi_gradient_correction_3d.f90 index 327422210..6773ad90b 100644 --- a/unit-tests/3d/test_mpi_gradient_correction_3d.f90 +++ b/unit-tests/3d/test_mpi_gradient_correction_3d.f90 @@ -118,7 +118,9 @@ program test_mpi_gradient_correction_3d parcels%B(4, :) = parcels%B(1, :) ! b33 needs to be calculated here - parcels%B(6, :) = get_b33(parcels%B(0:5, :), parcels%volume(:)) + do i =0, n_parcels + parcels%B(6, i) = get_b33(parcels%B(0:5, i), parcels%volume(i)) + enddo call vol2grid From 70b76edf6ef9e5a4a0543263f08c5e514d16efe9 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Rui=20Ap=C3=B3stolo?= Date: Wed, 8 Nov 2023 14:54:01 +0000 Subject: [PATCH 19/62] 1-indexed --- unit-tests/3d/test_mpi_gradient_correction_3d.f90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/unit-tests/3d/test_mpi_gradient_correction_3d.f90 b/unit-tests/3d/test_mpi_gradient_correction_3d.f90 index 6773ad90b..2d9f4ccc6 100644 --- a/unit-tests/3d/test_mpi_gradient_correction_3d.f90 +++ b/unit-tests/3d/test_mpi_gradient_correction_3d.f90 @@ -118,8 +118,8 @@ program test_mpi_gradient_correction_3d parcels%B(4, :) = parcels%B(1, :) ! b33 needs to be calculated here - do i =0, n_parcels - parcels%B(6, i) = get_b33(parcels%B(0:5, i), parcels%volume(i)) + do n = 1, n_parcels + parcels%B(6, n) = get_b33(parcels%B(0:5, n), parcels%volume(n)) enddo call vol2grid From 1d63ba729e9e2f9a2a4b52b4df51983ee7459123 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Rui=20Ap=C3=B3stolo?= Date: Wed, 8 Nov 2023 14:58:21 +0000 Subject: [PATCH 20/62] 1-indexed again --- unit-tests/3d/test_mpi_gradient_correction_3d.f90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/unit-tests/3d/test_mpi_gradient_correction_3d.f90 b/unit-tests/3d/test_mpi_gradient_correction_3d.f90 index 2d9f4ccc6..ab425b3e2 100644 --- a/unit-tests/3d/test_mpi_gradient_correction_3d.f90 +++ b/unit-tests/3d/test_mpi_gradient_correction_3d.f90 @@ -119,7 +119,7 @@ program test_mpi_gradient_correction_3d ! b33 needs to be calculated here do n = 1, n_parcels - parcels%B(6, n) = get_b33(parcels%B(0:5, n), parcels%volume(n)) + parcels%B(6, n) = get_b33(parcels%B(1:5, n), parcels%volume(n)) enddo call vol2grid From e3d93c519c9debcccfa453536861d6e262084019 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Rui=20Ap=C3=B3stolo?= Date: Mon, 13 Nov 2023 08:55:51 +0000 Subject: [PATCH 21/62] Initialised B33 in more tests --- unit-tests/3d/test_mpi_gradient_correction_3d.f90 | 2 +- unit-tests/3d/test_mpi_grid2par.f90 | 5 +++++ unit-tests/3d/test_mpi_parcel_diagnostics.f90 | 4 ++++ unit-tests/3d/test_mpi_parcel_init_3d.f90 | 1 + unit-tests/3d/test_mpi_parcel_pack.f90 | 8 +++++--- unit-tests/3d/test_mpi_parcel_split.f90 | 4 +++- unit-tests/3d/test_mpi_trilinear.f90 | 5 +++++ 7 files changed, 24 insertions(+), 5 deletions(-) diff --git a/unit-tests/3d/test_mpi_gradient_correction_3d.f90 b/unit-tests/3d/test_mpi_gradient_correction_3d.f90 index ab425b3e2..7be156473 100644 --- a/unit-tests/3d/test_mpi_gradient_correction_3d.f90 +++ b/unit-tests/3d/test_mpi_gradient_correction_3d.f90 @@ -117,7 +117,7 @@ program test_mpi_gradient_correction_3d ! b22 parcels%B(4, :) = parcels%B(1, :) - ! b33 needs to be calculated here + ! b33 do n = 1, n_parcels parcels%B(6, n) = get_b33(parcels%B(1:5, n), parcels%volume(n)) enddo diff --git a/unit-tests/3d/test_mpi_grid2par.f90 b/unit-tests/3d/test_mpi_grid2par.f90 index 6f72c1b6f..2141c968a 100644 --- a/unit-tests/3d/test_mpi_grid2par.f90 +++ b/unit-tests/3d/test_mpi_grid2par.f90 @@ -76,6 +76,11 @@ program test_mpi_grid2par ! b22 parcels%B(4, 1:n_parcels) = parcels%B(1, 1:n_parcels) + ! b33 + do n = 1, n_parcels + parcels%B(6, n) = get_b33(parcels%B(1:5, n), parcels%volume(n)) + enddo + velog(:, :, :, 1) = one velog(:, :, :, 2) = two velog(:, :, :, 3) = three diff --git a/unit-tests/3d/test_mpi_parcel_diagnostics.f90 b/unit-tests/3d/test_mpi_parcel_diagnostics.f90 index 3496b9a05..66db8e178 100644 --- a/unit-tests/3d/test_mpi_parcel_diagnostics.f90 +++ b/unit-tests/3d/test_mpi_parcel_diagnostics.f90 @@ -66,6 +66,10 @@ program test_mpi_parcel_diagnostics parcels%B(:, 1:n_parcels) = zero parcels%B(1, 1:n_parcels) = get_abc(parcels%volume(1:n_parcels)) ** f23 parcels%B(4, 1:n_parcels) = parcels%B(1, 1:n_parcels) + ! b33 + do n = 1, n_parcels + parcels%B(6, n) = get_b33(parcels%B(1:5, n), parcels%volume(n)) + enddo parcels%vorticity(:, 1:n_parcels) = f12 parcels%delta_pos(:, 1:n_parcels) = f12 diff --git a/unit-tests/3d/test_mpi_parcel_init_3d.f90 b/unit-tests/3d/test_mpi_parcel_init_3d.f90 index f3bc0effc..e0ae1d36a 100644 --- a/unit-tests/3d/test_mpi_parcel_init_3d.f90 +++ b/unit-tests/3d/test_mpi_parcel_init_3d.f90 @@ -108,6 +108,7 @@ program test_mpi_parcel_init_3d parcels%B(:, l) = zero parcels%B(1, l) = get_abc(v0) ** f23 parcels%B(4, l) = parcels%B(1, l) + parcels%B(6, l) = get_b33(parcels%B(1:5, l), parcels%volume(l)) l = l + 1 enddo enddo diff --git a/unit-tests/3d/test_mpi_parcel_pack.f90 b/unit-tests/3d/test_mpi_parcel_pack.f90 index 57dbb52bf..5800be3dd 100644 --- a/unit-tests/3d/test_mpi_parcel_pack.f90 +++ b/unit-tests/3d/test_mpi_parcel_pack.f90 @@ -47,10 +47,11 @@ program test_mpi_parcel_pack parcels%B(3, n) = 9.0d0 + a parcels%B(4, n) = 10.0d0 + a parcels%B(5, n) = 11.0d0 + a - parcels%volume(n) = 12.0d0 + a - parcels%buoyancy(n) = 13.0d0 + a + parcels%B(5, n) = 12.0d0 + a + parcels%volume(n) = 13.0d0 + a + parcels%buoyancy(n) = 14.0d0 + a #ifndef ENABLE_DRY_MODE - parcels%humidity(n) = 14.0d0 + a + parcels%humidity(n) = 15.0d0 + a #endif enddo @@ -92,6 +93,7 @@ program test_mpi_parcel_pack passed = (passed .and. (parcels%B(3, l) - buffer(i + IDX_B13) == zero)) passed = (passed .and. (parcels%B(4, l) - buffer(i + IDX_B22) == zero)) passed = (passed .and. (parcels%B(5, l) - buffer(i + IDX_B23) == zero)) + passed = (passed .and. (parcels%B(6, l) - buffer(i + IDX_B33) == zero)) passed = (passed .and. (parcels%volume(l) - buffer(i + IDX_VOL) == zero)) passed = (passed .and. (parcels%buoyancy(l) - buffer(i + IDX_BUO) == zero)) #ifndef ENABLE_DRY_MODE diff --git a/unit-tests/3d/test_mpi_parcel_split.f90 b/unit-tests/3d/test_mpi_parcel_split.f90 index e8766ec39..c75c648d3 100644 --- a/unit-tests/3d/test_mpi_parcel_split.f90 +++ b/unit-tests/3d/test_mpi_parcel_split.f90 @@ -141,7 +141,7 @@ program test_mpi_parcel_split subroutine fill_shape(m, angle) integer, intent(in) :: m double precision, intent(in) :: angle - double precision :: B11, B12, B13, B22, B23, st, ct, sp, cp + double precision :: B11, B12, B13, B22, B23, B33, st, ct, sp, cp st = dsin(angle) ct = dcos(angle) @@ -153,12 +153,14 @@ subroutine fill_shape(m, angle) B13 = (a2 - c2) * ct * sp * cp B22 = a2 * st ** 2 * sp ** 2 + b2 * ct ** 2 + c2 * st ** 2 * cp ** 2 B23 = (a2 - c2) * st * sp * cp + B33 = a2 * cp ** 2 + c2 * sp ** 2 parcels%B(1, m) = B11 parcels%B(2, m) = B12 parcels%B(3, m) = B13 parcels%B(4, m) = B22 parcels%B(5, m) = B23 + parcels%B(6, m) = B33 end subroutine fill_shape end program test_mpi_parcel_split diff --git a/unit-tests/3d/test_mpi_trilinear.f90 b/unit-tests/3d/test_mpi_trilinear.f90 index 28636d908..9c3a3109e 100644 --- a/unit-tests/3d/test_mpi_trilinear.f90 +++ b/unit-tests/3d/test_mpi_trilinear.f90 @@ -78,6 +78,11 @@ program test_mpi_trilinear ! b22 parcels%B(4, 1:n_parcels) = parcels%B(1, 1:n_parcels) + ! b33 + do n = 1, n_parcels + parcels%B(6, n) = get_b33(parcels%B(1:5, n), parcels%volume(n)) + enddo + call par2grid error = abs(get_sum(volg) - dble(ngrid) * vcell) From 49c9bebde6e6eac9d44007dc6858d771ff55f2d2 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Rui=20Ap=C3=B3stolo?= Date: Mon, 13 Nov 2023 09:00:26 +0000 Subject: [PATCH 22/62] Updated imports on test_mpi_parcel_init_3d.f90 --- unit-tests/3d/test_mpi_parcel_init_3d.f90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/unit-tests/3d/test_mpi_parcel_init_3d.f90 b/unit-tests/3d/test_mpi_parcel_init_3d.f90 index e0ae1d36a..69d8a5aba 100644 --- a/unit-tests/3d/test_mpi_parcel_init_3d.f90 +++ b/unit-tests/3d/test_mpi_parcel_init_3d.f90 @@ -12,7 +12,7 @@ program test_mpi_parcel_init_3d use parcel_container use parcel_init, only : init_timer, parcel_default, init_parcels_from_grids use parcel_interpl, only : par2grid, par2grid_timer, halo_swap_timer - use parcel_ellipsoid, only : get_abc + use parcel_ellipsoid, only : get_abc, get_b33 use fields, only : tbuoyg, field_default use field_ops, only : get_rms, get_abs_max use parameters, only : update_parameters, dx, nx, ny, nz, lower, vcell From 85362253596925f8f186e3a10ec03a02c0c30a32 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Rui=20Ap=C3=B3stolo?= Date: Mon, 13 Nov 2023 09:02:36 +0000 Subject: [PATCH 23/62] Updated variable declaration on test_mpi_trilinear.f90 --- unit-tests/3d/test_mpi_trilinear.f90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/unit-tests/3d/test_mpi_trilinear.f90 b/unit-tests/3d/test_mpi_trilinear.f90 index 9c3a3109e..f9b7af3fb 100644 --- a/unit-tests/3d/test_mpi_trilinear.f90 +++ b/unit-tests/3d/test_mpi_trilinear.f90 @@ -18,7 +18,7 @@ program test_mpi_trilinear implicit none double precision :: error - integer :: ix, iy, iz, i, j, k, l, n_per_dim + integer :: ix, iy, iz, i, j, k, l, n_per_dim, n double precision :: im, corner(3) logical :: passed = .true. From e2315559934e33d2fe5c50163cf902a5b686a6eb Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Rui=20Ap=C3=B3stolo?= Date: Mon, 13 Nov 2023 09:04:27 +0000 Subject: [PATCH 24/62] Updated imports on test_mpi_trilinear.f90 --- unit-tests/3d/test_mpi_trilinear.f90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/unit-tests/3d/test_mpi_trilinear.f90 b/unit-tests/3d/test_mpi_trilinear.f90 index f9b7af3fb..92dc568ef 100644 --- a/unit-tests/3d/test_mpi_trilinear.f90 +++ b/unit-tests/3d/test_mpi_trilinear.f90 @@ -10,7 +10,7 @@ program test_mpi_trilinear use parcel_container use mpi_layout use parcel_interpl, only : par2grid, par2grid_timer, halo_swap_timer - use parcel_ellipsoid, only : get_abc + use parcel_ellipsoid, only : get_abc, get_b33 use parameters, only : lower, update_parameters, vcell, dx, nx, ny, nz, ngrid use fields, only : volg, field_alloc use field_ops, only : get_sum From 0e17cb4963bed1e644855a49ff554b5d01a6663f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Rui=20Ap=C3=B3stolo?= Date: Mon, 13 Nov 2023 09:06:22 +0000 Subject: [PATCH 25/62] Updated import and variable declaration on test_mpi_grid2par.f90 --- unit-tests/3d/test_mpi_grid2par.f90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/unit-tests/3d/test_mpi_grid2par.f90 b/unit-tests/3d/test_mpi_grid2par.f90 index 2141c968a..4b2d7f873 100644 --- a/unit-tests/3d/test_mpi_grid2par.f90 +++ b/unit-tests/3d/test_mpi_grid2par.f90 @@ -11,14 +11,14 @@ program test_mpi_grid2par use constants, only : pi, zero, one, two, three, four, five, f12, f23 use parcel_container use parcel_interpl, only : grid2par, grid2par_timer - use parcel_ellipsoid, only : get_abc + use parcel_ellipsoid, only : get_abc, get_b33 use parameters, only : lower, update_parameters, vcell, dx, nx, ny, nz use fields, only : velog, vortg, velgradg, field_alloc use mpi_timer implicit none double precision :: error - integer :: ix, iy, iz, i, j, k, l, n_per_dim + integer :: ix, iy, iz, i, j, k, l, n_per_dim, n double precision :: im, corner(3) logical :: passed = .true. From 64bd0a9f2c8899945ba1161abbde30a60df8d0d2 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Rui=20Ap=C3=B3stolo?= Date: Mon, 13 Nov 2023 09:12:45 +0000 Subject: [PATCH 26/62] Updated variable declaration on test_mpi_parcel_diagnostic.f90 --- unit-tests/3d/test_mpi_parcel_diagnostics.f90 | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/unit-tests/3d/test_mpi_parcel_diagnostics.f90 b/unit-tests/3d/test_mpi_parcel_diagnostics.f90 index 66db8e178..9d0c07f72 100644 --- a/unit-tests/3d/test_mpi_parcel_diagnostics.f90 +++ b/unit-tests/3d/test_mpi_parcel_diagnostics.f90 @@ -8,13 +8,14 @@ program test_mpi_parcel_diagnostics use mpi_layout use parcel_container use parcel_diagnostics + use parcel_ellipsoid, only: get_b33 use parameters, only : lower, update_parameters, extent, nx, ny, nz, vcell, dx, set_vmin use mpi_timer implicit none logical :: passed = .true. integer, parameter :: n_per_dim = 2 - integer :: ix, iy, iz, i, j, k, l, n_total + integer :: ix, iy, iz, i, j, k, l, n_total, n double precision :: im, corner(3), total_vol call mpi_env_initialise From cbb290d406168ab68beb16e9be58614a291142d6 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Rui=20Ap=C3=B3stolo?= Date: Mon, 13 Nov 2023 09:27:58 +0000 Subject: [PATCH 27/62] Updated test_mpi_parcel_correction_3d.f90 and test_mpi_parecl_pack.f90 --- unit-tests/3d/test_mpi_parcel_correction_3d.f90 | 7 ++++++- unit-tests/3d/test_mpi_parcel_pack.f90 | 2 +- 2 files changed, 7 insertions(+), 2 deletions(-) diff --git a/unit-tests/3d/test_mpi_parcel_correction_3d.f90 b/unit-tests/3d/test_mpi_parcel_correction_3d.f90 index 01a391251..950ed9289 100644 --- a/unit-tests/3d/test_mpi_parcel_correction_3d.f90 +++ b/unit-tests/3d/test_mpi_parcel_correction_3d.f90 @@ -18,7 +18,7 @@ program test_parcel_correction_3d , vort_corr_timer & , init_parcel_correction use parcel_interpl, only : vol2grid, halo_swap_timer - use parcel_ellipsoid, only : get_abc + use parcel_ellipsoid, only : get_abc, get_b33 use parcel_init, only : init_regular_positions use parameters, only : lower, extent, update_parameters, vcell, nx, ny, nz, dx use fields, only : volg, field_default @@ -120,6 +120,11 @@ program test_parcel_correction_3d ! b22 parcels%B(4, :) = parcels%B(1, :) + ! b33 + do n = 1, n_parcels + parcels%B(6, n) = get_b33(parcels%B(1:5, n), parcels%volume(n)) + enddo + call vol2grid volg(0:nz, lo(2):hi(2), lo(1):hi(1)) = abs(volg(0:nz, lo(2):hi(2), lo(1):hi(1)) / vcell - one) diff --git a/unit-tests/3d/test_mpi_parcel_pack.f90 b/unit-tests/3d/test_mpi_parcel_pack.f90 index 5800be3dd..bcfedd62e 100644 --- a/unit-tests/3d/test_mpi_parcel_pack.f90 +++ b/unit-tests/3d/test_mpi_parcel_pack.f90 @@ -47,7 +47,7 @@ program test_mpi_parcel_pack parcels%B(3, n) = 9.0d0 + a parcels%B(4, n) = 10.0d0 + a parcels%B(5, n) = 11.0d0 + a - parcels%B(5, n) = 12.0d0 + a + parcels%B(6, n) = 12.0d0 + a parcels%volume(n) = 13.0d0 + a parcels%buoyancy(n) = 14.0d0 + a #ifndef ENABLE_DRY_MODE From a8fc1b3824a055571b4fa6a94fb42b45e7b275e8 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Rui=20Ap=C3=B3stolo?= Date: Mon, 13 Nov 2023 10:05:43 +0000 Subject: [PATCH 28/62] get_dBdt needs Bout(IDX_B33) fixed --- src/3d/stepper/rk_utils.f90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/3d/stepper/rk_utils.f90 b/src/3d/stepper/rk_utils.f90 index 0fb6551f1..211afad95 100644 --- a/src/3d/stepper/rk_utils.f90 +++ b/src/3d/stepper/rk_utils.f90 @@ -75,7 +75,7 @@ function get_dBdt(Bin, S, vorticity, volume) result(Bout) - S(I_DUDX) * Bin(I_B23) & ! - du/dx * B23 + dvdz * Bin(I_B33) ! + dv/dz * B33 - Bout(I_B33) = get_B33(Bout(I_B11:I_B23), volume) + Bout(I_B33) = end function get_dBdt ! Estimate a suitable time step based on the velocity strain From 76739f9d74acf93d7b3e76bccf9b8d465835d386 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Rui=20Ap=C3=B3stolo?= Date: Mon, 13 Nov 2023 11:24:56 +0000 Subject: [PATCH 29/62] Changed get_dBdt --- src/3d/stepper/rk_utils.f90 | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/src/3d/stepper/rk_utils.f90 b/src/3d/stepper/rk_utils.f90 index 211afad95..941357268 100644 --- a/src/3d/stepper/rk_utils.f90 +++ b/src/3d/stepper/rk_utils.f90 @@ -75,7 +75,10 @@ function get_dBdt(Bin, S, vorticity, volume) result(Bout) - S(I_DUDX) * Bin(I_B23) & ! - du/dx * B23 + dvdz * Bin(I_B33) ! + dv/dz * B33 - Bout(I_B33) = + + ! dB33/dt + Bout(I_B33) = two * (S(I_DWDX) * Bin(I_B13) + S(I_DWDY) * Bin(I_B23) + dwdz * Bin(I_B33)) + end function get_dBdt ! Estimate a suitable time step based on the velocity strain From f9aeb9053906d8ee6096765a2cef3b7a20cd904f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Rui=20Ap=C3=B3stolo?= Date: Mon, 13 Nov 2023 11:28:36 +0000 Subject: [PATCH 30/62] Changed get_dBdt to remove volume --- src/3d/stepper/ls_rk.f90 | 6 ++---- src/3d/stepper/rk_utils.f90 | 3 +-- 2 files changed, 3 insertions(+), 6 deletions(-) diff --git a/src/3d/stepper/ls_rk.f90 b/src/3d/stepper/ls_rk.f90 index 5e73c0e66..03493c13d 100644 --- a/src/3d/stepper/ls_rk.f90 +++ b/src/3d/stepper/ls_rk.f90 @@ -143,8 +143,7 @@ subroutine ls_rk_substep(dt, step) do n = 1, n_parcels parcels%delta_b(:, n) = get_dBdt(parcels%B(:, n), & parcels%strain(:, n), & - parcels%vorticity(:, n), & - parcels%volume(n)) + parcels%vorticity(:, n)) enddo !$omp end parallel do @@ -163,8 +162,7 @@ subroutine ls_rk_substep(dt, step) parcels%delta_b(:, n) = parcels%delta_b(:, n) & + get_dBdt(parcels%B(:, n), & parcels%strain(:, n), & - parcels%vorticity(:, n), & - parcels%volume(n)) + parcels%vorticity(:, n)) enddo !$omp end parallel do diff --git a/src/3d/stepper/rk_utils.f90 b/src/3d/stepper/rk_utils.f90 index 941357268..070cb1c0f 100644 --- a/src/3d/stepper/rk_utils.f90 +++ b/src/3d/stepper/rk_utils.f90 @@ -26,11 +26,10 @@ module rk_utils ! @param[in] vorticity of parcel ! @param[in] volume is the parcel volume ! @returns dB/dt in Bout - function get_dBdt(Bin, S, vorticity, volume) result(Bout) + function get_dBdt(Bin, S, vorticity) result(Bout) double precision, intent(in) :: Bin(I_B33) double precision, intent(in) :: S(5) double precision, intent(in) :: vorticity(n_dim) - double precision, intent(in) :: volume double precision :: Bout(I_B33) double precision :: dudz, dvdx, dvdz, dwdz From 4c8c029e18b4db469fe7285db2423fa9ccb50089 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Rui=20Ap=C3=B3stolo?= Date: Mon, 13 Nov 2023 12:09:26 +0000 Subject: [PATCH 31/62] test debug --- unit-tests/3d/test_mpi_parcel_pack.f90 | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/unit-tests/3d/test_mpi_parcel_pack.f90 b/unit-tests/3d/test_mpi_parcel_pack.f90 index bcfedd62e..72f6f36c7 100644 --- a/unit-tests/3d/test_mpi_parcel_pack.f90 +++ b/unit-tests/3d/test_mpi_parcel_pack.f90 @@ -34,6 +34,7 @@ program test_mpi_parcel_pack call parcel_alloc(n_parcels) ! initialise parcels + write (*, *) "initialise parcels" do n = 1, n_parcels a = dble(n-1) * 100.0d0 parcels%position(1, n) = 1.0d0 + a @@ -55,6 +56,8 @@ program test_mpi_parcel_pack #endif enddo + write (*, *) "init done" + do m = 1, 100 ! randomly select parcels to pack pid = -1 @@ -70,12 +73,15 @@ program test_mpi_parcel_pack pid(k) = l enddo + write (*, *) "allocate(buffer(n_par_attrib * n_pack))" allocate(buffer(n_par_attrib * n_pack)) ! pack the parcels + write (*, *) "call parcel_pack(pid, n_pack, buffer)" call parcel_pack(pid, n_pack, buffer) ! check the buffer of packed parcels + write (*, *) "check the buffer of packed parcels" do n = 1, n_pack l = pid(n) @@ -104,6 +110,7 @@ program test_mpi_parcel_pack deallocate(buffer) enddo + write (*, *) " mpi checks" if (world%rank == world%root) then call MPI_Reduce(MPI_IN_PLACE, passed, 1, MPI_LOGICAL, MPI_LAND, world%root, world%comm, world%err) else From b7163f28aba703fa0a49ad239d458f3a60012ced Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Rui=20Ap=C3=B3stolo?= Date: Mon, 13 Nov 2023 13:21:36 +0000 Subject: [PATCH 32/62] Revert "test debug" This reverts commit 4c8c029e18b4db469fe7285db2423fa9ccb50089. --- unit-tests/3d/test_mpi_parcel_pack.f90 | 7 ------- 1 file changed, 7 deletions(-) diff --git a/unit-tests/3d/test_mpi_parcel_pack.f90 b/unit-tests/3d/test_mpi_parcel_pack.f90 index 72f6f36c7..bcfedd62e 100644 --- a/unit-tests/3d/test_mpi_parcel_pack.f90 +++ b/unit-tests/3d/test_mpi_parcel_pack.f90 @@ -34,7 +34,6 @@ program test_mpi_parcel_pack call parcel_alloc(n_parcels) ! initialise parcels - write (*, *) "initialise parcels" do n = 1, n_parcels a = dble(n-1) * 100.0d0 parcels%position(1, n) = 1.0d0 + a @@ -56,8 +55,6 @@ program test_mpi_parcel_pack #endif enddo - write (*, *) "init done" - do m = 1, 100 ! randomly select parcels to pack pid = -1 @@ -73,15 +70,12 @@ program test_mpi_parcel_pack pid(k) = l enddo - write (*, *) "allocate(buffer(n_par_attrib * n_pack))" allocate(buffer(n_par_attrib * n_pack)) ! pack the parcels - write (*, *) "call parcel_pack(pid, n_pack, buffer)" call parcel_pack(pid, n_pack, buffer) ! check the buffer of packed parcels - write (*, *) "check the buffer of packed parcels" do n = 1, n_pack l = pid(n) @@ -110,7 +104,6 @@ program test_mpi_parcel_pack deallocate(buffer) enddo - write (*, *) " mpi checks" if (world%rank == world%root) then call MPI_Reduce(MPI_IN_PLACE, passed, 1, MPI_LOGICAL, MPI_LAND, world%root, world%comm, world%err) else From 45e8f009f9c2ed4d896c793b3e1df69cefd63b16 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Rui=20Ap=C3=B3stolo?= Date: Mon, 13 Nov 2023 13:23:30 +0000 Subject: [PATCH 33/62] Fixed test_mpi_laplace_correction --- unit-tests/3d/test_mpi_laplace_correction_3d.f90 | 7 ++++++- 1 file changed, 6 insertions(+), 1 deletion(-) diff --git a/unit-tests/3d/test_mpi_laplace_correction_3d.f90 b/unit-tests/3d/test_mpi_laplace_correction_3d.f90 index ed508427a..4d3257025 100644 --- a/unit-tests/3d/test_mpi_laplace_correction_3d.f90 +++ b/unit-tests/3d/test_mpi_laplace_correction_3d.f90 @@ -16,7 +16,7 @@ program test_mpi_laplace_correction_3d , vort_corr_timer & , init_parcel_correction use parcel_interpl, only : vol2grid, halo_swap_timer - use parcel_ellipsoid, only : get_abc + use parcel_ellipsoid, only : get_abc, get_b33 use parcel_init, only : init_regular_positions use parameters, only : lower, extent, update_parameters, vcell, nx, ny, nz, dx use fields, only : volg, field_default @@ -117,6 +117,11 @@ program test_mpi_laplace_correction_3d ! b22 parcels%B(4, :) = parcels%B(1, :) + ! b33 + do n = 1, n_parcels + parcels%B(6, n) = get_b33(parcels%B(1:5, n), parcels%volume(n)) + enddo + call vol2grid volg(0:nz, lo(2):hi(2), lo(1):hi(1)) = abs(volg(0:nz, lo(2):hi(2), lo(1):hi(1)) / vcell - one) From 063803a3b58949c2a63b4f956b2d23aa3e298d0b Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Rui=20Ap=C3=B3stolo?= Date: Mon, 27 Nov 2023 10:13:34 +0000 Subject: [PATCH 34/62] Formula for get_dBdt Bout(I_B33) --- src/3d/stepper/rk_utils.f90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/3d/stepper/rk_utils.f90 b/src/3d/stepper/rk_utils.f90 index 070cb1c0f..f62ef7d0a 100644 --- a/src/3d/stepper/rk_utils.f90 +++ b/src/3d/stepper/rk_utils.f90 @@ -75,7 +75,7 @@ function get_dBdt(Bin, S, vorticity) result(Bout) + dvdz * Bin(I_B33) ! + dv/dz * B33 - ! dB33/dt + ! dB33/dt = 2 * (dw/dx * B13 + dw/dy * B23 + dw/dz * B33) Bout(I_B33) = two * (S(I_DWDX) * Bin(I_B13) + S(I_DWDY) * Bin(I_B23) + dwdz * Bin(I_B33)) end function get_dBdt From 3c454f8c536613d53fd6808e57e34d10d5749f41 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Rui=20Ap=C3=B3stolo?= Date: Wed, 6 Dec 2023 14:24:26 +0000 Subject: [PATCH 35/62] Fixed Bm(6:n) in parcel_merge.f90 --- src/3d/parcels/parcel_merge.f90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/3d/parcels/parcel_merge.f90 b/src/3d/parcels/parcel_merge.f90 index 3820130e6..f1d776c50 100644 --- a/src/3d/parcels/parcel_merge.f90 +++ b/src/3d/parcels/parcel_merge.f90 @@ -271,7 +271,7 @@ subroutine do_group_merge(isma, iclo, n_merge, Bm, vm) Bm(3, n) = Bm(3, n) + mu * (five * delx * delz + parcels%B(3, is)) Bm(4, n) = Bm(4, n) + mu * (five * dely ** 2 + parcels%B(4, is)) Bm(5, n) = Bm(5, n) + mu * (five * dely * delz + parcels%B(5, is)) - Bm(6, n) = Bm(6, n) + mu * (five * delz ** 2 + parcels%B(5, is)) + Bm(6, n) = Bm(6, n) + mu * (five * delz ** 2 + parcels%B(6, is)) enddo end subroutine do_group_merge From f352f14fb2949464a76ece92c0a47964774d1806 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Rui=20Ap=C3=B3stolo?= Date: Fri, 12 Jan 2024 14:04:28 +0000 Subject: [PATCH 36/62] B matrix normalisation at parcel initialisation --- src/3d/parcels/parcel_merge.f90 | 1 - src/3d/stepper/ls_rk.f90 | 14 ++++++++++++++ 2 files changed, 14 insertions(+), 1 deletion(-) diff --git a/src/3d/parcels/parcel_merge.f90 b/src/3d/parcels/parcel_merge.f90 index f1d776c50..97d2c4753 100644 --- a/src/3d/parcels/parcel_merge.f90 +++ b/src/3d/parcels/parcel_merge.f90 @@ -13,7 +13,6 @@ module parcel_merging , get_dely_across_periodic & , parcel_delete use parcel_ellipsoid, only : get_abc - ! use parcel_ellipsoid, only : get_B33, get_abc use options, only : parcel #if defined (ENABLE_VERBOSE) && !defined (NDEBUG) use options, only : verbose diff --git a/src/3d/stepper/ls_rk.f90 b/src/3d/stepper/ls_rk.f90 index 03493c13d..140ed1d29 100644 --- a/src/3d/stepper/ls_rk.f90 +++ b/src/3d/stepper/ls_rk.f90 @@ -3,11 +3,13 @@ ! (see https://doi.org/10.5194/gmd-10-3145-2017) ! ============================================================================= module ls_rk + use constants, only : f13 use options, only : time use dimensions, only : I_Z use parcel_container use parcel_bc use parcel_mpi, only : parcel_communicate + use parcel_ellipsoid, only : get_abc use rk_utils, only: get_dBdt, get_time_step use utils, only : write_step use parcel_interpl, only : par2grid, grid2par @@ -131,6 +133,7 @@ subroutine ls_rk_substep(dt, step) double precision, intent(in) :: dt integer, intent(in) :: step double precision :: ca, cb + double precision :: factor, detB integer :: n ca = captr(step) @@ -163,9 +166,11 @@ subroutine ls_rk_substep(dt, step) + get_dBdt(parcels%B(:, n), & parcels%strain(:, n), & parcels%vorticity(:, n)) + enddo !$omp end parallel do + call stop_timer(rk_timer) endif @@ -180,6 +185,15 @@ subroutine ls_rk_substep(dt, step) + cb * dt * parcels%delta_vor(:, n) parcels%B(:, n) = parcels%B(:, n) & + cb * dt * parcels%delta_b(:, n) + + ! normalize B matrix + detB = parcels%B(1, n) * (parcels%B(4, n) * parcels%B(6, n) - parcels%B(5, n) ** 2) & + - parcels%B(2, n) * (parcels%B(2, n) * parcels%B(6, n) - parcels%B(3, n) * parcels%B(5, n)) & + + parcels%B(3, n) * (parcels%B(2, n) * parcels%B(5, n) - parcels%B(3, n) * parcels%B(4, n)) + + factor = (get_abc(parcels%volume(n)) ** 2 / detB) ** f13 + + parcels%B(:, n) = parcels%B(1:6, n) * factor enddo !$omp end parallel do From ef7f685912f1cfb717a51915dacd4e12858c9d2f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Rui=20Ap=C3=B3stolo?= Date: Mon, 15 Jan 2024 10:09:51 +0000 Subject: [PATCH 37/62] B33 init in parcel_init.f90 --- src/3d/parcels/parcel_init.f90 | 3 +++ 1 file changed, 3 insertions(+) diff --git a/src/3d/parcels/parcel_init.f90 b/src/3d/parcels/parcel_init.f90 index ae66cf325..7602aaf75 100644 --- a/src/3d/parcels/parcel_init.f90 +++ b/src/3d/parcels/parcel_init.f90 @@ -89,6 +89,9 @@ subroutine parcel_default ! B22 parcels%B(4, n) = l23 + + ! B33 + parcels%B(6, n) = l23 enddo !$omp end do !$omp end parallel From dd57f70136542a9810e1dbcb6ab6cea93b9f4767 Mon Sep 17 00:00:00 2001 From: sjboeing Date: Thu, 18 Jan 2024 15:44:51 +0000 Subject: [PATCH 38/62] Parcel merge and split bugfixes --- src/3d/parcels/parcel_merge.f90 | 2 +- src/3d/parcels/parcel_split.f90 | 1 + 2 files changed, 2 insertions(+), 1 deletion(-) diff --git a/src/3d/parcels/parcel_merge.f90 b/src/3d/parcels/parcel_merge.f90 index 97d2c4753..21e6b6d27 100644 --- a/src/3d/parcels/parcel_merge.f90 +++ b/src/3d/parcels/parcel_merge.f90 @@ -312,7 +312,7 @@ subroutine geometric_merge(isma, iclo, n_merge) factor = (get_abc(V(l)) ** 2 / detB) ** f13 - parcels%B(:, ic) = B(1:5, l) * factor + parcels%B(:, ic) = B(:, l) * factor endif enddo diff --git a/src/3d/parcels/parcel_split.f90 b/src/3d/parcels/parcel_split.f90 index b2f9ba178..8fc2f443f 100644 --- a/src/3d/parcels/parcel_split.f90 +++ b/src/3d/parcels/parcel_split.f90 @@ -129,6 +129,7 @@ subroutine parcel_split parcels%B(3, n) = B(3) - f34 * D(1) * V(1, 1) * V(3, 1) parcels%B(4, n) = B(4) - f34 * D(1) * V(2, 1) ** 2 parcels%B(5, n) = B(5) - f34 * D(1) * V(2, 1) * V(3, 1) + parcels%B(6, n) = B(6) - f34 * D(1) * V(3, 1) ** 2 parcels%volume(n) = f12 * vol From 04004cad408d8dfc952b446a3400fe1bf4fa7401 Mon Sep 17 00:00:00 2001 From: sjboeing Date: Fri, 19 Jan 2024 16:52:54 +0000 Subject: [PATCH 39/62] RK module: small fixes to (mostly) OpenMP --- src/3d/stepper/ls_rk.f90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/3d/stepper/ls_rk.f90 b/src/3d/stepper/ls_rk.f90 index 140ed1d29..1c9f87af8 100644 --- a/src/3d/stepper/ls_rk.f90 +++ b/src/3d/stepper/ls_rk.f90 @@ -176,7 +176,7 @@ subroutine ls_rk_substep(dt, step) call start_timer(rk_timer) - !$omp parallel do default(shared) private(n) + !$omp parallel do default(shared) private(n,detB,factor) do n = 1, n_parcels parcels%position(:, n) = parcels%position(:, n) & + cb * dt * parcels%delta_pos(:, n) @@ -193,7 +193,7 @@ subroutine ls_rk_substep(dt, step) factor = (get_abc(parcels%volume(n)) ** 2 / detB) ** f13 - parcels%B(:, n) = parcels%B(1:6, n) * factor + parcels%B(:, n) = parcels%B(:, n) * factor enddo !$omp end parallel do From 6c663e1b776a37c2438f66753969e38b40213a79 Mon Sep 17 00:00:00 2001 From: sjboeing Date: Sun, 21 Jan 2024 09:31:28 +0000 Subject: [PATCH 40/62] Correct delta_b to be consistent with actual B increment --- src/3d/stepper/ls_rk.f90 | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/src/3d/stepper/ls_rk.f90 b/src/3d/stepper/ls_rk.f90 index 1c9f87af8..934f2d54d 100644 --- a/src/3d/stepper/ls_rk.f90 +++ b/src/3d/stepper/ls_rk.f90 @@ -3,7 +3,7 @@ ! (see https://doi.org/10.5194/gmd-10-3145-2017) ! ============================================================================= module ls_rk - use constants, only : f13 + use constants, only : f13, one use options, only : time use dimensions, only : I_Z use parcel_container @@ -183,6 +183,7 @@ subroutine ls_rk_substep(dt, step) parcels%vorticity(:, n) = parcels%vorticity(:, n) & + cb * dt * parcels%delta_vor(:, n) + parcels%B(:, n) = parcels%B(:, n) & + cb * dt * parcels%delta_b(:, n) @@ -193,6 +194,7 @@ subroutine ls_rk_substep(dt, step) factor = (get_abc(parcels%volume(n)) ** 2 / detB) ** f13 + parcels%delta_b(:, n) = parcels%delta_b(:, n) + parcels%B(:, n) * (factor - one)/(cb * dt) parcels%B(:, n) = parcels%B(:, n) * factor enddo !$omp end parallel do From 3b6d994ef0d4f6dc2a7554359602136b5600244f Mon Sep 17 00:00:00 2001 From: sjboeing Date: Mon, 22 Jan 2024 11:43:12 +0000 Subject: [PATCH 41/62] Small alignment correction --- src/3d/stepper/ls_rk.f90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/3d/stepper/ls_rk.f90 b/src/3d/stepper/ls_rk.f90 index 934f2d54d..d30c67c1e 100644 --- a/src/3d/stepper/ls_rk.f90 +++ b/src/3d/stepper/ls_rk.f90 @@ -203,7 +203,7 @@ subroutine ls_rk_substep(dt, step) if (step == n_stages) then call parcel_communicate - return + return endif call start_timer(rk_timer) From 36e6b0d1f6e40529160865e77d1d173525303d61 Mon Sep 17 00:00:00 2001 From: sjboeing Date: Sun, 4 Feb 2024 17:39:00 +0000 Subject: [PATCH 42/62] Use gridded fields for dbdt --- src/3d/fields/fields.f90 | 3 +++ src/3d/parcels/parcel_container.f90 | 16 +++++++++----- src/3d/parcels/parcel_interpl.f90 | 24 ++++++++++++++++---- src/3d/stepper/ls_rk.f90 | 6 ++--- src/3d/stepper/rk_utils.f90 | 34 ++++++++++------------------- 5 files changed, 48 insertions(+), 35 deletions(-) diff --git a/src/3d/fields/fields.f90 b/src/3d/fields/fields.f90 index 2550fbf90..17e51a269 100644 --- a/src/3d/fields/fields.f90 +++ b/src/3d/fields/fields.f90 @@ -52,7 +52,10 @@ module fields ! velocity strain indices integer, parameter :: I_DUDX = 1 & ! index for du/dx strain component , I_DUDY = 2 & ! index for du/dy strain component + , I_DUDZ = 6 & ! index for du/dz strain component + , I_DVDX = 7 & ! index for dv/dy strain component , I_DVDY = 3 & ! index for dv/dy strain component + , I_DVDZ = 8 & ! index for dv/dy strain component , I_DWDX = 4 & ! index for dw/dx strain component , I_DWDY = 5 ! index for dw/dy strain component diff --git a/src/3d/parcels/parcel_container.f90 b/src/3d/parcels/parcel_container.f90 index 2302f10ab..dd905e7fc 100644 --- a/src/3d/parcels/parcel_container.f90 +++ b/src/3d/parcels/parcel_container.f90 @@ -56,7 +56,10 @@ module parcel_container IDX_RK4_DB33, & ! RK4 variable for B33 IDX_RK4_DUDX, & ! RK4 variable du/dx IDX_RK4_DUDY, & ! RK4 variable du/dy + IDX_RK4_DUDZ, & ! RK4 variable du/dz + IDX_RK4_DVDX, & ! RK4 variable dv/dx IDX_RK4_DVDY, & ! RK4 variable dv/dy + IDX_RK4_DVDZ, & ! RK4 variable dv/dz IDX_RK4_DWDX, & ! RK4 variable dw/dx IDX_RK4_DWDY ! RK4 variable dw/dy @@ -134,11 +137,14 @@ subroutine set_buffer_indices IDX_RK4_DB33 = i + 11 IDX_RK4_DUDX = i + 12 IDX_RK4_DUDY = i + 13 - IDX_RK4_DVDY = i + 14 - IDX_RK4_DWDX = i + 15 - IDX_RK4_DWDY = i + 16 + IDX_RK4_DUDZ = i + 14 + IDX_RK4_DVDX = i + 15 + IDX_RK4_DVDY = i + 16 + IDX_RK4_DVDZ = i + 17 + IDX_RK4_DWDX = i + 18 + IDX_RK4_DWDY = i + 19 - i = i + 17 + i = i + 20 n_par_attrib = set_ellipsoid_buffer_indices(i) @@ -319,7 +325,7 @@ subroutine parcel_alloc(num) ! LS-RK4 variables allocate(parcels%delta_pos(3, num)) allocate(parcels%delta_vor(3, num)) - allocate(parcels%strain(5, num)) + allocate(parcels%strain(8, num)) allocate(parcels%delta_b(6, num)) end subroutine parcel_alloc diff --git a/src/3d/parcels/parcel_interpl.f90 b/src/3d/parcels/parcel_interpl.f90 index 0b9082297..852c7be3e 100644 --- a/src/3d/parcels/parcel_interpl.f90 +++ b/src/3d/parcels/parcel_interpl.f90 @@ -433,10 +433,26 @@ subroutine grid2par(add) parcels%delta_pos(l, n) = parcels%delta_pos(l, n) & + f14 * sum(weights * velog(ks:ks+1, js:js+1, is:is+1, l)) enddo - do l = 1,5 - parcels%strain(l, n) = parcels%strain(l, n) & - + f14 * sum(weights * velgradg(ks:ks+1, js:js+1, is:is+1, l)) - enddo + parcels%strain(I_DUDX, n) = parcels%strain(I_DUDX, n) & + + f14 * sum(weights * velgradg(ks:ks+1, js:js+1, is:is+1, I_DUDX)) + parcels%strain(I_DUDY, n) = parcels%strain(I_DUDY, n) & + + f14 * sum(weights * velgradg(ks:ks+1, js:js+1, is:is+1, I_DUDY)) + parcels%strain(I_DUDZ, n) = parcels%strain(I_DUDZ, n) & + + f14 * sum(weights * (vortg(ks:ks+1, js:js+1, is:is+1, I_Y) + & + velgradg(ks:ks+1, js:js+1, is:is+1, I_DWDX))) + parcels%strain(I_DVDX, n) = parcels%strain(I_DVDX, n) & + + f14 * sum(weights * (vortg(ks:ks+1, js:js+1, is:is+1, I_Z) + & + velgradg(ks:ks+1, js:js+1, is:is+1, I_DUDY))) + parcels%strain(I_DVDY, n) = parcels%strain(I_DVDY, n) & + + f14 * sum(weights * velgradg(ks:ks+1, js:js+1, is:is+1, I_DVDY)) + parcels%strain(I_DVDZ, n) = parcels%strain(I_DVDZ, n) & + + f14 * sum(weights * ( -vortg(ks:ks+1, js:js+1, is:is+1,I_X) + & + velgradg(ks:ks+1, js:js+1, is:is+1, I_DWDY))) + parcels%strain(I_DWDX, n) = parcels%strain(I_DWDX, n) & + + f14 * sum(weights * velgradg(ks:ks+1, js:js+1, is:is+1, I_DWDX)) + parcels%strain(I_DWDY, n) = parcels%strain(I_DWDY, n) & + + f14 * sum(weights * velgradg(ks:ks+1, js:js+1, is:is+1, I_DWDY)) + do l = 1,3 parcels%delta_vor(l, n) = parcels%delta_vor(l, n) & + f14 * sum(weights * vtend(ks:ks+1, js:js+1, is:is+1, l)) diff --git a/src/3d/stepper/ls_rk.f90 b/src/3d/stepper/ls_rk.f90 index d30c67c1e..52fda8481 100644 --- a/src/3d/stepper/ls_rk.f90 +++ b/src/3d/stepper/ls_rk.f90 @@ -145,8 +145,7 @@ subroutine ls_rk_substep(dt, step) !$omp parallel do default(shared) private(n) do n = 1, n_parcels parcels%delta_b(:, n) = get_dBdt(parcels%B(:, n), & - parcels%strain(:, n), & - parcels%vorticity(:, n)) + parcels%strain(:, n)) enddo !$omp end parallel do @@ -164,8 +163,7 @@ subroutine ls_rk_substep(dt, step) do n = 1, n_parcels parcels%delta_b(:, n) = parcels%delta_b(:, n) & + get_dBdt(parcels%B(:, n), & - parcels%strain(:, n), & - parcels%vorticity(:, n)) + parcels%strain(:, n)) enddo !$omp end parallel do diff --git a/src/3d/stepper/rk_utils.f90 b/src/3d/stepper/rk_utils.f90 index f62ef7d0a..00426d502 100644 --- a/src/3d/stepper/rk_utils.f90 +++ b/src/3d/stepper/rk_utils.f90 @@ -1,7 +1,7 @@ module rk_utils use dimensions, only : n_dim, I_X, I_Y, I_Z use parcel_ellipsoid, only : get_B33, I_B11, I_B12, I_B13, I_B22, I_B23, I_B33 - use fields, only : velgradg, tbuoyg, vortg, I_DUDX, I_DUDY, I_DVDY, I_DWDX, I_DWDY + use fields, only : velgradg, tbuoyg, vortg, I_DUDX, I_DUDY, I_DUDZ, I_DVDX, I_DVDY, I_DVDZ, I_DWDX, I_DWDY use field_mpi, only : field_halo_fill use constants, only : zero, one, two, f12 use parameters, only : nx, ny, nz, dxi, vcell @@ -26,21 +26,11 @@ module rk_utils ! @param[in] vorticity of parcel ! @param[in] volume is the parcel volume ! @returns dB/dt in Bout - function get_dBdt(Bin, S, vorticity) result(Bout) + function get_dBdt(Bin, S) result(Bout) double precision, intent(in) :: Bin(I_B33) - double precision, intent(in) :: S(5) - double precision, intent(in) :: vorticity(n_dim) + double precision, intent(in) :: S(8) double precision :: Bout(I_B33) - double precision :: dudz, dvdx, dvdz, dwdz - - ! du/dz = \eta + dw/dx - dudz = vorticity(I_Y) + S(I_DWDX) - - ! dv/dx \zeta + du/dy - dvdx = vorticity(I_Z) + S(I_DUDY) - - ! dv/dz = dw/dy - \xi - dvdz = S(I_DWDY) - vorticity(I_X) + double precision :: dwdz ! dw/dz = - (du/dx + dv/dy) dwdz = - (S(I_DUDX) + S(I_DVDY)) @@ -48,31 +38,31 @@ function get_dBdt(Bin, S, vorticity) result(Bout) ! B33 = get_B33(Bin, volume) ! dB11/dt = 2 * (du/dx * B11 + du/dy * B12 + du/dz * B13) - Bout(I_B11) = two * (S(I_DUDX) * Bin(I_B11) + S(I_DUDY) * Bin(I_B12) + dudz * Bin(I_B13)) + Bout(I_B11) = two * (S(I_DUDX) * Bin(I_B11) + S(I_DUDY) * Bin(I_B12) + S(I_DUDZ) * Bin(I_B13)) ! dB12/dt = - Bout(I_B12) = dvdx * Bin(I_B11) & ! dv/dx * B11 + Bout(I_B12) = S(I_DVDX) * Bin(I_B11) & ! dv/dx * B11 - dwdz * Bin(I_B12) & ! - dw/dz * B12 - + dvdz * Bin(I_B13) & ! + dv/dz * B13 + + S(I_DVDZ) * Bin(I_B13) & ! + dv/dz * B13 + S(I_DUDY) * Bin(I_B22) & ! + du/dy * B22 - + dudz * Bin(I_B23) ! + du/dz * B23 + + S(I_DUDZ) * Bin(I_B23) ! + du/dz * B23 ! dB13/dt = Bout(I_B13) = S(I_DWDX) * Bin(I_B11) & ! dw/dx * B11 + S(I_DWDY) * Bin(I_B12) & ! + dw/dy * B12 - S(I_DVDY) * Bin(I_B13) & ! - dv/dy * B13 + S(I_DUDY) * Bin(I_B23) & ! + du/dy * B23 - + dudz * Bin(I_B33) ! + du/dz * B33 + + S(I_DUDZ) * Bin(I_B33) ! + du/dz * B33 ! dB22/dt = 2 * (dv/dx * B12 + dv/dy * B22 + dv/dz * B23) - Bout(I_B22) = two * (dvdx * Bin(I_B12) + S(I_DVDY) * Bin(I_B22) + dvdz * Bin(I_B23)) + Bout(I_B22) = two * (S(I_DVDX) * Bin(I_B12) + S(I_DVDY) * Bin(I_B22) + S(I_DVDZ) * Bin(I_B23)) ! dB23/dt = Bout(I_B23) = S(I_DWDX) * Bin(I_B12) & ! dw/dx * B12 - + dvdx * Bin(I_B13) & ! + dv/dx * B13 + + S(I_DVDX) * Bin(I_B13) & ! + dv/dx * B13 + S(I_DWDY) * Bin(I_B22) & ! + dw/dy * B22 - S(I_DUDX) * Bin(I_B23) & ! - du/dx * B23 - + dvdz * Bin(I_B33) ! + dv/dz * B33 + + S(I_DVDZ) * Bin(I_B33) ! + dv/dz * B33 ! dB33/dt = 2 * (dw/dx * B13 + dw/dy * B23 + dw/dz * B33) From b9b013931d92e39916e558c8ca54b946b6165cda Mon Sep 17 00:00:00 2001 From: Steven Boeing Date: Mon, 5 Feb 2024 16:30:23 +0000 Subject: [PATCH 43/62] Work on matrix exponential branch --- src/3d/parcels/parcel_interpl.f90 | 6 +- src/3d/stepper/ls_rk.f90 | 25 +------- src/3d/stepper/rk_utils.f90 | 101 +++++++++++++++--------------- 3 files changed, 58 insertions(+), 74 deletions(-) diff --git a/src/3d/parcels/parcel_interpl.f90 b/src/3d/parcels/parcel_interpl.f90 index 852c7be3e..35b2e78c1 100644 --- a/src/3d/parcels/parcel_interpl.f90 +++ b/src/3d/parcels/parcel_interpl.f90 @@ -400,6 +400,7 @@ subroutine grid2par(add) do n = 1, n_parcels parcels%delta_pos(:, n) = zero parcels%delta_vor(:, n) = zero + parcels%strain(:, n) = zero enddo !$omp end do !$omp end parallel @@ -410,6 +411,7 @@ subroutine grid2par(add) do n = 1, n_parcels parcels%delta_pos(:, n) = zero parcels%delta_vor(:, n) = zero + parcels%strain(:, n) = zero enddo !$omp end do !$omp end parallel @@ -419,8 +421,6 @@ subroutine grid2par(add) !$omp do private(n, l, p, points, is, js, ks, weights) do n = 1, n_parcels - parcels%strain(:, n) = zero - points = get_ellipsoid_points(parcels%position(:, n), & parcels%B(:, n), n) @@ -452,7 +452,7 @@ subroutine grid2par(add) + f14 * sum(weights * velgradg(ks:ks+1, js:js+1, is:is+1, I_DWDX)) parcels%strain(I_DWDY, n) = parcels%strain(I_DWDY, n) & + f14 * sum(weights * velgradg(ks:ks+1, js:js+1, is:is+1, I_DWDY)) - + do l = 1,3 parcels%delta_vor(l, n) = parcels%delta_vor(l, n) & + f14 * sum(weights * vtend(ks:ks+1, js:js+1, is:is+1, l)) diff --git a/src/3d/stepper/ls_rk.f90 b/src/3d/stepper/ls_rk.f90 index 52fda8481..091679bb5 100644 --- a/src/3d/stepper/ls_rk.f90 +++ b/src/3d/stepper/ls_rk.f90 @@ -10,7 +10,7 @@ module ls_rk use parcel_bc use parcel_mpi, only : parcel_communicate use parcel_ellipsoid, only : get_abc - use rk_utils, only: get_dBdt, get_time_step + use rk_utils, only: get_time_step, evolve_ellipsoid use utils, only : write_step use parcel_interpl, only : par2grid, grid2par use fields, only : velgradg, velog, vortg, vtend, tbuoyg @@ -142,13 +142,6 @@ subroutine ls_rk_substep(dt, step) if (step == 1) then call start_timer(rk_timer) - !$omp parallel do default(shared) private(n) - do n = 1, n_parcels - parcels%delta_b(:, n) = get_dBdt(parcels%B(:, n), & - parcels%strain(:, n)) - enddo - !$omp end parallel do - call stop_timer(rk_timer) else call vor2vel @@ -159,16 +152,6 @@ subroutine ls_rk_substep(dt, step) call start_timer(rk_timer) - !$omp parallel do default(shared) private(n) - do n = 1, n_parcels - parcels%delta_b(:, n) = parcels%delta_b(:, n) & - + get_dBdt(parcels%B(:, n), & - parcels%strain(:, n)) - - enddo - !$omp end parallel do - - call stop_timer(rk_timer) endif @@ -182,8 +165,7 @@ subroutine ls_rk_substep(dt, step) parcels%vorticity(:, n) = parcels%vorticity(:, n) & + cb * dt * parcels%delta_vor(:, n) - parcels%B(:, n) = parcels%B(:, n) & - + cb * dt * parcels%delta_b(:, n) + evolve_ellipsoid(parcels%B(:, n), parcels%strain(:, n), cb * dt) ! normalize B matrix detB = parcels%B(1, n) * (parcels%B(4, n) * parcels%B(6, n) - parcels%B(5, n) ** 2) & @@ -192,7 +174,6 @@ subroutine ls_rk_substep(dt, step) factor = (get_abc(parcels%volume(n)) ** 2 / detB) ** f13 - parcels%delta_b(:, n) = parcels%delta_b(:, n) + parcels%B(:, n) * (factor - one)/(cb * dt) parcels%B(:, n) = parcels%B(:, n) * factor enddo !$omp end parallel do @@ -210,7 +191,7 @@ subroutine ls_rk_substep(dt, step) do n = 1, n_parcels parcels%delta_pos(:, n) = ca * parcels%delta_pos(:, n) parcels%delta_vor(:, n) = ca * parcels%delta_vor(:, n) - parcels%delta_b(:, n) = ca * parcels%delta_b(:, n) + parcels%strain(:, n) = ca * parcels%stain(:, n) enddo !$omp end parallel do diff --git a/src/3d/stepper/rk_utils.f90 b/src/3d/stepper/rk_utils.f90 index 00426d502..b250a087a 100644 --- a/src/3d/stepper/rk_utils.f90 +++ b/src/3d/stepper/rk_utils.f90 @@ -20,55 +20,58 @@ module rk_utils contains - ! Advance the B matrix. - ! @param[in] Bin are the B matrix components of the parcel - ! @param[in] S is the local velocity strain - ! @param[in] vorticity of parcel - ! @param[in] volume is the parcel volume - ! @returns dB/dt in Bout - function get_dBdt(Bin, S) result(Bout) - double precision, intent(in) :: Bin(I_B33) - double precision, intent(in) :: S(8) - double precision :: Bout(I_B33) - double precision :: dwdz - - ! dw/dz = - (du/dx + dv/dy) - dwdz = - (S(I_DUDX) + S(I_DVDY)) - - ! B33 = get_B33(Bin, volume) - - ! dB11/dt = 2 * (du/dx * B11 + du/dy * B12 + du/dz * B13) - Bout(I_B11) = two * (S(I_DUDX) * Bin(I_B11) + S(I_DUDY) * Bin(I_B12) + S(I_DUDZ) * Bin(I_B13)) - - ! dB12/dt = - Bout(I_B12) = S(I_DVDX) * Bin(I_B11) & ! dv/dx * B11 - - dwdz * Bin(I_B12) & ! - dw/dz * B12 - + S(I_DVDZ) * Bin(I_B13) & ! + dv/dz * B13 - + S(I_DUDY) * Bin(I_B22) & ! + du/dy * B22 - + S(I_DUDZ) * Bin(I_B23) ! + du/dz * B23 - - ! dB13/dt = - Bout(I_B13) = S(I_DWDX) * Bin(I_B11) & ! dw/dx * B11 - + S(I_DWDY) * Bin(I_B12) & ! + dw/dy * B12 - - S(I_DVDY) * Bin(I_B13) & ! - dv/dy * B13 - + S(I_DUDY) * Bin(I_B23) & ! + du/dy * B23 - + S(I_DUDZ) * Bin(I_B33) ! + du/dz * B33 - - ! dB22/dt = 2 * (dv/dx * B12 + dv/dy * B22 + dv/dz * B23) - Bout(I_B22) = two * (S(I_DVDX) * Bin(I_B12) + S(I_DVDY) * Bin(I_B22) + S(I_DVDZ) * Bin(I_B23)) - - ! dB23/dt = - Bout(I_B23) = S(I_DWDX) * Bin(I_B12) & ! dw/dx * B12 - + S(I_DVDX) * Bin(I_B13) & ! + dv/dx * B13 - + S(I_DWDY) * Bin(I_B22) & ! + dw/dy * B22 - - S(I_DUDX) * Bin(I_B23) & ! - du/dx * B23 - + S(I_DVDZ) * Bin(I_B33) ! + dv/dz * B33 - - - ! dB33/dt = 2 * (dw/dx * B13 + dw/dy * B23 + dw/dz * B33) - Bout(I_B33) = two * (S(I_DWDX) * Bin(I_B13) + S(I_DWDY) * Bin(I_B23) + dwdz * Bin(I_B33)) - - end function get_dBdt + subroutine evolve_ellipsoid(B, S, dt_sub) +- double precision, intent(inout) :: Bin(I_B33) +- double precision, intent(in) :: S(8) + double precision, intent(in) :: dt_sub + double precision :: Bmat(3,3) + double precision :: Smat(3,3) + double precision :: Qmat(3,3) + double precision :: Imat(3,3) + + Bmat(1, 1) = B(I_B11) ! B11 + Bmat(1, 2) = B(I_B12) ! B12 + Bmat(1, 3) = B(I_B13) ! B13 + Bmat(2, 1) = B(I_B12) ! B21 + Bmat(2, 2) = B(I_B22) ! B22 + Bmat(2, 3) = B(I_B23) ! B23 + Bmat(3, 1) = B(I_B13) ! B31 + Bmat(3, 2) = B(I_B23) ! B32 + Bmat(3, 3) = B(I_B33) ! B33 + + Smat(1, 1) = S(I_DUDX) ! S11 + Smat(1, 2) = S(I_DUDY) ! S12 + Smat(1, 3) = S(I_DUDZ) ! S13 + Smat(2, 1) = S(I_DVDX) ! S21 + Smat(2, 2) = S(I_DVDY) ! S22 + Smat(2, 3) = S(I_DVDZ) ! S23 + Smat(3, 1) = S(I_DWDX) ! S31 + Smat(3, 2) = S(I_DWDY) ! S32 + Smat(3, 3) = -(S(I_DUDX) + S(I_DVDY)) ! S33 + + Imat = 0.0 + Imat(1, 1) = 1.0 + Imat(2, 2) = 1.0 + Imat(3, 3) = 1.0 + + Qmat = Imat + & + 0.25 * dt_sub * Smat + & + 0.03125 * dt_sub * dt_sub * matmul(Smat, Smat) + & + (1.0 / 384.0) * dt_sub * dt_sub * dt_sub * matmul(Smat, matmul(Smat, Smat)) + + Qmat = matmul(Qmat, Qmat) + Qmat = matmul(Qmat, Qmat) + Bmat = matmul(Qmat, matmul(Bmat, transpose(Qmat))) + + B(I_B11) = Bmat(1, 1) + B(I_B12) = Bmat(1, 2) + B(I_B13) = Bmat(1, 3) + B(I_B22) = Bmat(2, 2) + B(I_B23) = Bmat(2, 3) + B(I_B33) = Bmat(3, 3) + + end subroutine + ! Estimate a suitable time step based on the velocity strain ! and buoyancy gradient. From 9c2309b6fa1cf9a6f406d6ef5d26f6af4a78092d Mon Sep 17 00:00:00 2001 From: sjboeing Date: Mon, 5 Feb 2024 20:23:58 +0000 Subject: [PATCH 44/62] Fix bugs, remove delta_b --- src/3d/parcels/parcel_container.f90 | 36 ++++++++++------------------- src/3d/stepper/ls_rk.f90 | 8 +++---- src/3d/stepper/rk_utils.f90 | 11 +++------ 3 files changed, 19 insertions(+), 36 deletions(-) diff --git a/src/3d/parcels/parcel_container.f90 b/src/3d/parcels/parcel_container.f90 index dd905e7fc..3fa18c3f2 100644 --- a/src/3d/parcels/parcel_container.f90 +++ b/src/3d/parcels/parcel_container.f90 @@ -86,8 +86,8 @@ module parcel_container double precision, allocatable, dimension(:, :) :: & delta_pos, & ! velocity delta_vor, & ! vorticity tendency - strain, & - delta_b ! B-matrix tendency + strain + end type parcel_container_type type(parcel_container_type) parcels @@ -129,22 +129,16 @@ subroutine set_buffer_indices IDX_RK4_X_DVOR = i + 3 IDX_RK4_Y_DVOR = i + 4 IDX_RK4_Z_DVOR = i + 5 - IDX_RK4_DB11 = i + 6 - IDX_RK4_DB12 = i + 7 - IDX_RK4_DB13 = i + 8 - IDX_RK4_DB22 = i + 9 - IDX_RK4_DB23 = i + 10 - IDX_RK4_DB33 = i + 11 - IDX_RK4_DUDX = i + 12 - IDX_RK4_DUDY = i + 13 - IDX_RK4_DUDZ = i + 14 - IDX_RK4_DVDX = i + 15 - IDX_RK4_DVDY = i + 16 - IDX_RK4_DVDZ = i + 17 - IDX_RK4_DWDX = i + 18 - IDX_RK4_DWDY = i + 19 - - i = i + 20 + IDX_RK4_DUDX = i + 6 + IDX_RK4_DUDY = i + 7 + IDX_RK4_DUDZ = i + 8 + IDX_RK4_DVDX = i + 9 + IDX_RK4_DVDY = i + 10 + IDX_RK4_DVDZ = i + 11 + IDX_RK4_DWDX = i + 12 + IDX_RK4_DWDY = i + 13 + + i = i + 14 n_par_attrib = set_ellipsoid_buffer_indices(i) @@ -261,7 +255,6 @@ subroutine parcel_replace(n, m) ! LS-RK4 variables: parcels%delta_pos(:, n) = parcels%delta_pos(:, m) parcels%delta_vor(:, n) = parcels%delta_vor(:, m) - parcels%delta_b(:, n) = parcels%delta_b(:, m) parcels%strain(:, n) = parcels%strain(:, m) end subroutine parcel_replace @@ -297,7 +290,6 @@ subroutine parcel_resize(new_size) call resize_array(parcels%delta_pos, new_size, n_parcels) call resize_array(parcels%delta_vor, new_size, n_parcels) call resize_array(parcels%strain, new_size, n_parcels) - call resize_array(parcels%delta_b, new_size, n_parcels) call stop_timer(resize_timer) @@ -326,7 +318,6 @@ subroutine parcel_alloc(num) allocate(parcels%delta_pos(3, num)) allocate(parcels%delta_vor(3, num)) allocate(parcels%strain(8, num)) - allocate(parcels%delta_b(6, num)) end subroutine parcel_alloc @@ -356,7 +347,6 @@ subroutine parcel_dealloc deallocate(parcels%delta_pos) deallocate(parcels%delta_vor) deallocate(parcels%strain) - deallocate(parcels%delta_b) end subroutine parcel_dealloc @@ -378,7 +368,6 @@ subroutine parcel_serialize(n, buffer) ! LS-RK4 variables: buffer(IDX_RK4_X_DPOS:IDX_RK4_Z_DPOS) = parcels%delta_pos(:, n) buffer(IDX_RK4_X_DVOR:IDX_RK4_Z_DVOR) = parcels%delta_vor(:, n) - buffer(IDX_RK4_DB11:IDX_RK4_DB33) = parcels%delta_b(:, n) buffer(IDX_RK4_DUDX:IDX_RK4_DWDY) = parcels%strain(:, n) call parcel_ellipsoid_serialize(n, buffer) @@ -402,7 +391,6 @@ subroutine parcel_deserialize(n, buffer) ! LS-RK4 variables: parcels%delta_pos(:, n) = buffer(IDX_RK4_X_DPOS:IDX_RK4_Z_DPOS) parcels%delta_vor(:, n) = buffer(IDX_RK4_X_DVOR:IDX_RK4_Z_DVOR) - parcels%delta_b(:, n) = buffer(IDX_RK4_DB11:IDX_RK4_DB33) parcels%strain(:, n) = buffer(IDX_RK4_DUDX:IDX_RK4_DWDY) call parcel_ellipsoid_deserialize(n, buffer) diff --git a/src/3d/stepper/ls_rk.f90 b/src/3d/stepper/ls_rk.f90 index 091679bb5..ad032dc07 100644 --- a/src/3d/stepper/ls_rk.f90 +++ b/src/3d/stepper/ls_rk.f90 @@ -114,6 +114,7 @@ subroutine ls_rk_step(t) enddo call ls_rk_substep(dt, n_stages) + ! normalize B matrix call start_timer(rk_timer) call apply_parcel_reflective_bc call stop_timer(rk_timer) @@ -157,7 +158,7 @@ subroutine ls_rk_substep(dt, step) call start_timer(rk_timer) - !$omp parallel do default(shared) private(n,detB,factor) + !$omp parallel do default(shared) private(n, detB, factor) do n = 1, n_parcels parcels%position(:, n) = parcels%position(:, n) & + cb * dt * parcels%delta_pos(:, n) @@ -165,9 +166,8 @@ subroutine ls_rk_substep(dt, step) parcels%vorticity(:, n) = parcels%vorticity(:, n) & + cb * dt * parcels%delta_vor(:, n) - evolve_ellipsoid(parcels%B(:, n), parcels%strain(:, n), cb * dt) + call evolve_ellipsoid(parcels%B(:, n), parcels%strain(:, n), cb * dt) - ! normalize B matrix detB = parcels%B(1, n) * (parcels%B(4, n) * parcels%B(6, n) - parcels%B(5, n) ** 2) & - parcels%B(2, n) * (parcels%B(2, n) * parcels%B(6, n) - parcels%B(3, n) * parcels%B(5, n)) & + parcels%B(3, n) * (parcels%B(2, n) * parcels%B(5, n) - parcels%B(3, n) * parcels%B(4, n)) @@ -191,7 +191,7 @@ subroutine ls_rk_substep(dt, step) do n = 1, n_parcels parcels%delta_pos(:, n) = ca * parcels%delta_pos(:, n) parcels%delta_vor(:, n) = ca * parcels%delta_vor(:, n) - parcels%strain(:, n) = ca * parcels%stain(:, n) + parcels%strain(:, n) = ca * parcels%strain(:, n) enddo !$omp end parallel do diff --git a/src/3d/stepper/rk_utils.f90 b/src/3d/stepper/rk_utils.f90 index b250a087a..548acf418 100644 --- a/src/3d/stepper/rk_utils.f90 +++ b/src/3d/stepper/rk_utils.f90 @@ -17,17 +17,17 @@ module rk_utils #endif implicit none + double precision, parameter :: Imat(3,3)=reshape((/one,zero,zero,zero,one,zero,zero,zero,one/), (/3,3/)) contains subroutine evolve_ellipsoid(B, S, dt_sub) -- double precision, intent(inout) :: Bin(I_B33) -- double precision, intent(in) :: S(8) + double precision, intent(inout) :: B(I_B33) + double precision, intent(in) :: S(8) double precision, intent(in) :: dt_sub double precision :: Bmat(3,3) double precision :: Smat(3,3) double precision :: Qmat(3,3) - double precision :: Imat(3,3) Bmat(1, 1) = B(I_B11) ! B11 Bmat(1, 2) = B(I_B12) ! B12 @@ -49,11 +49,6 @@ subroutine evolve_ellipsoid(B, S, dt_sub) Smat(3, 2) = S(I_DWDY) ! S32 Smat(3, 3) = -(S(I_DUDX) + S(I_DVDY)) ! S33 - Imat = 0.0 - Imat(1, 1) = 1.0 - Imat(2, 2) = 1.0 - Imat(3, 3) = 1.0 - Qmat = Imat + & 0.25 * dt_sub * Smat + & 0.03125 * dt_sub * dt_sub * matmul(Smat, Smat) + & From f178974a8b1a1e4d659d84668ba58dd5a3ff0282 Mon Sep 17 00:00:00 2001 From: sjboeing Date: Tue, 6 Feb 2024 10:16:30 +0000 Subject: [PATCH 45/62] Using one more iteration in the exponentiation --- src/3d/stepper/rk_utils.f90 | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/src/3d/stepper/rk_utils.f90 b/src/3d/stepper/rk_utils.f90 index 548acf418..d1fd63506 100644 --- a/src/3d/stepper/rk_utils.f90 +++ b/src/3d/stepper/rk_utils.f90 @@ -50,10 +50,11 @@ subroutine evolve_ellipsoid(B, S, dt_sub) Smat(3, 3) = -(S(I_DUDX) + S(I_DVDY)) ! S33 Qmat = Imat + & - 0.25 * dt_sub * Smat + & - 0.03125 * dt_sub * dt_sub * matmul(Smat, Smat) + & - (1.0 / 384.0) * dt_sub * dt_sub * dt_sub * matmul(Smat, matmul(Smat, Smat)) + 0.125 * dt_sub * Smat + & + 0.0078125 * dt_sub * dt_sub * matmul(Smat, Smat) + & + (one / 3072.0) * dt_sub * dt_sub * dt_sub * matmul(Smat, matmul(Smat, Smat)) + Qmat = matmul(Qmat, Qmat) Qmat = matmul(Qmat, Qmat) Qmat = matmul(Qmat, Qmat) Bmat = matmul(Qmat, matmul(Bmat, transpose(Qmat))) From a841d813666ad9a8e95b135cff09327f7f5bbf8a Mon Sep 17 00:00:00 2001 From: sjboeing Date: Tue, 6 Feb 2024 10:19:27 +0000 Subject: [PATCH 46/62] Using velgradg with 8 components --- src/3d/fields/fields.f90 | 14 +++++------ src/3d/inversion/inversion.f90 | 8 +++++++ src/3d/parcels/parcel_container.f90 | 30 ++++++++++++------------ src/3d/parcels/parcel_interpl.f90 | 36 ++++++++--------------------- src/3d/stepper/ls_rk.f90 | 4 ++-- 5 files changed, 42 insertions(+), 50 deletions(-) diff --git a/src/3d/fields/fields.f90 b/src/3d/fields/fields.f90 index 17e51a269..2f43f92d9 100644 --- a/src/3d/fields/fields.f90 +++ b/src/3d/fields/fields.f90 @@ -52,12 +52,12 @@ module fields ! velocity strain indices integer, parameter :: I_DUDX = 1 & ! index for du/dx strain component , I_DUDY = 2 & ! index for du/dy strain component - , I_DUDZ = 6 & ! index for du/dz strain component - , I_DVDX = 7 & ! index for dv/dy strain component - , I_DVDY = 3 & ! index for dv/dy strain component - , I_DVDZ = 8 & ! index for dv/dy strain component - , I_DWDX = 4 & ! index for dw/dx strain component - , I_DWDY = 5 ! index for dw/dy strain component + , I_DUDZ = 3 & ! index for du/dz strain component + , I_DVDX = 4 & ! index for dv/dy strain component + , I_DVDY = 5 & ! index for dv/dy strain component + , I_DVDZ = 6 & ! index for dv/dy strain component + , I_DWDX = 7 & ! index for dw/dx strain component + , I_DWDY = 8 ! index for dw/dy strain component contains @@ -77,7 +77,7 @@ subroutine field_alloc hhi = box%hhi allocate(velog(hlo(3):hhi(3), hlo(2):hhi(2), hlo(1):hhi(1), n_dim)) - allocate(velgradg(hlo(3):hhi(3), hlo(2):hhi(2), hlo(1):hhi(1), 5)) + allocate(velgradg(hlo(3):hhi(3), hlo(2):hhi(2), hlo(1):hhi(1), 8)) allocate(volg(hlo(3):hhi(3), hlo(2):hhi(2), hlo(1):hhi(1))) diff --git a/src/3d/inversion/inversion.f90 b/src/3d/inversion/inversion.f90 index e56a9940e..c7c4e968f 100644 --- a/src/3d/inversion/inversion.f90 +++ b/src/3d/inversion/inversion.f90 @@ -307,6 +307,14 @@ subroutine vel2vgrad(svel) call field_halo_fill(velgradg, l_alloc=.true.) + !$omp parallel workshare + ! fill the other components + velgradg(:, :, :, I_DUDZ) = velgradg(:, :, :, I_DWDX) + vortg(:, : , :, I_Y) + velgradg(:, :, :, I_DVDX) = velgradg(:, :, :, I_DUDY) + vortg(:, : , :, I_Z) + velgradg(:, :, :, I_DVDZ) = velgradg(:, :, :, I_DWDY) - vortg(:, : , :, I_X) + !$omp end parallel workshare + + end subroutine vel2vgrad !:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: diff --git a/src/3d/parcels/parcel_container.f90 b/src/3d/parcels/parcel_container.f90 index 3fa18c3f2..1e871fb13 100644 --- a/src/3d/parcels/parcel_container.f90 +++ b/src/3d/parcels/parcel_container.f90 @@ -86,7 +86,7 @@ module parcel_container double precision, allocatable, dimension(:, :) :: & delta_pos, & ! velocity delta_vor, & ! vorticity tendency - strain + int_strain end type parcel_container_type @@ -241,21 +241,21 @@ subroutine parcel_replace(n, m) endif #endif - parcels%position(:, n) = parcels%position(:, m) - parcels%vorticity(:, n) = parcels%vorticity(:, m) - parcels%volume(n) = parcels%volume(m) - parcels%buoyancy(n) = parcels%buoyancy(m) + parcels%position(:, n) = parcels%position(:, m) + parcels%vorticity(:, n) = parcels%vorticity(:, m) + parcels%volume(n) = parcels%volume(m) + parcels%buoyancy(n) = parcels%buoyancy(m) #ifndef ENABLE_DRY_MODE - parcels%humidity(n) = parcels%humidity(m) + parcels%humidity(n) = parcels%humidity(m) #endif - parcels%B(:, n) = parcels%B(:, m) + parcels%B(:, n) = parcels%B(:, m) call parcel_ellipsoid_replace(n, m) ! LS-RK4 variables: - parcels%delta_pos(:, n) = parcels%delta_pos(:, m) - parcels%delta_vor(:, n) = parcels%delta_vor(:, m) - parcels%strain(:, n) = parcels%strain(:, m) + parcels%delta_pos(:, n) = parcels%delta_pos(:, m) + parcels%delta_vor(:, n) = parcels%delta_vor(:, m) + parcels%int_strain(:, n) = parcels%int_strain(:, m) end subroutine parcel_replace @@ -289,7 +289,7 @@ subroutine parcel_resize(new_size) ! LS-RK4 variables call resize_array(parcels%delta_pos, new_size, n_parcels) call resize_array(parcels%delta_vor, new_size, n_parcels) - call resize_array(parcels%strain, new_size, n_parcels) + call resize_array(parcels%int_strain, new_size, n_parcels) call stop_timer(resize_timer) @@ -317,7 +317,7 @@ subroutine parcel_alloc(num) ! LS-RK4 variables allocate(parcels%delta_pos(3, num)) allocate(parcels%delta_vor(3, num)) - allocate(parcels%strain(8, num)) + allocate(parcels%int_strain(8, num)) end subroutine parcel_alloc @@ -346,7 +346,7 @@ subroutine parcel_dealloc ! LS-RK4 variables deallocate(parcels%delta_pos) deallocate(parcels%delta_vor) - deallocate(parcels%strain) + deallocate(parcels%int_strain) end subroutine parcel_dealloc @@ -368,7 +368,7 @@ subroutine parcel_serialize(n, buffer) ! LS-RK4 variables: buffer(IDX_RK4_X_DPOS:IDX_RK4_Z_DPOS) = parcels%delta_pos(:, n) buffer(IDX_RK4_X_DVOR:IDX_RK4_Z_DVOR) = parcels%delta_vor(:, n) - buffer(IDX_RK4_DUDX:IDX_RK4_DWDY) = parcels%strain(:, n) + buffer(IDX_RK4_DUDX:IDX_RK4_DWDY) = parcels%int_strain(:, n) call parcel_ellipsoid_serialize(n, buffer) end subroutine parcel_serialize @@ -391,7 +391,7 @@ subroutine parcel_deserialize(n, buffer) ! LS-RK4 variables: parcels%delta_pos(:, n) = buffer(IDX_RK4_X_DPOS:IDX_RK4_Z_DPOS) parcels%delta_vor(:, n) = buffer(IDX_RK4_X_DVOR:IDX_RK4_Z_DVOR) - parcels%strain(:, n) = buffer(IDX_RK4_DUDX:IDX_RK4_DWDY) + parcels%int_strain(:, n) = buffer(IDX_RK4_DUDX:IDX_RK4_DWDY) call parcel_ellipsoid_deserialize(n, buffer) end subroutine parcel_deserialize diff --git a/src/3d/parcels/parcel_interpl.f90 b/src/3d/parcels/parcel_interpl.f90 index 35b2e78c1..43d1d041f 100644 --- a/src/3d/parcels/parcel_interpl.f90 +++ b/src/3d/parcels/parcel_interpl.f90 @@ -400,7 +400,7 @@ subroutine grid2par(add) do n = 1, n_parcels parcels%delta_pos(:, n) = zero parcels%delta_vor(:, n) = zero - parcels%strain(:, n) = zero + parcels%int_strain(:, n) = zero enddo !$omp end do !$omp end parallel @@ -411,7 +411,7 @@ subroutine grid2par(add) do n = 1, n_parcels parcels%delta_pos(:, n) = zero parcels%delta_vor(:, n) = zero - parcels%strain(:, n) = zero + parcels%int_strain(:, n) = zero enddo !$omp end do !$omp end parallel @@ -430,32 +430,16 @@ subroutine grid2par(add) ! loop over grid points which are part of the interpolation do l = 1,3 - parcels%delta_pos(l, n) = parcels%delta_pos(l, n) & - + f14 * sum(weights * velog(ks:ks+1, js:js+1, is:is+1, l)) + parcels%delta_pos(l, n) = parcels%delta_pos(l, n) & + + f14 * sum(weights * velog(ks:ks+1, js:js+1, is:is+1, l)) + enddo + do l = 1, 8 + parcels%int_strain(l, n) = parcels%int_strain(l, n) & + + f14 * sum(weights * velgradg(ks:ks+1, js:js+1, is:is+1, l)) enddo - parcels%strain(I_DUDX, n) = parcels%strain(I_DUDX, n) & - + f14 * sum(weights * velgradg(ks:ks+1, js:js+1, is:is+1, I_DUDX)) - parcels%strain(I_DUDY, n) = parcels%strain(I_DUDY, n) & - + f14 * sum(weights * velgradg(ks:ks+1, js:js+1, is:is+1, I_DUDY)) - parcels%strain(I_DUDZ, n) = parcels%strain(I_DUDZ, n) & - + f14 * sum(weights * (vortg(ks:ks+1, js:js+1, is:is+1, I_Y) + & - velgradg(ks:ks+1, js:js+1, is:is+1, I_DWDX))) - parcels%strain(I_DVDX, n) = parcels%strain(I_DVDX, n) & - + f14 * sum(weights * (vortg(ks:ks+1, js:js+1, is:is+1, I_Z) + & - velgradg(ks:ks+1, js:js+1, is:is+1, I_DUDY))) - parcels%strain(I_DVDY, n) = parcels%strain(I_DVDY, n) & - + f14 * sum(weights * velgradg(ks:ks+1, js:js+1, is:is+1, I_DVDY)) - parcels%strain(I_DVDZ, n) = parcels%strain(I_DVDZ, n) & - + f14 * sum(weights * ( -vortg(ks:ks+1, js:js+1, is:is+1,I_X) + & - velgradg(ks:ks+1, js:js+1, is:is+1, I_DWDY))) - parcels%strain(I_DWDX, n) = parcels%strain(I_DWDX, n) & - + f14 * sum(weights * velgradg(ks:ks+1, js:js+1, is:is+1, I_DWDX)) - parcels%strain(I_DWDY, n) = parcels%strain(I_DWDY, n) & - + f14 * sum(weights * velgradg(ks:ks+1, js:js+1, is:is+1, I_DWDY)) - do l = 1,3 - parcels%delta_vor(l, n) = parcels%delta_vor(l, n) & - + f14 * sum(weights * vtend(ks:ks+1, js:js+1, is:is+1, l)) + parcels%delta_vor(l, n) = parcels%delta_vor(l, n) & + + f14 * sum(weights * vtend(ks:ks+1, js:js+1, is:is+1, l)) enddo enddo enddo diff --git a/src/3d/stepper/ls_rk.f90 b/src/3d/stepper/ls_rk.f90 index ad032dc07..068a0513d 100644 --- a/src/3d/stepper/ls_rk.f90 +++ b/src/3d/stepper/ls_rk.f90 @@ -166,7 +166,7 @@ subroutine ls_rk_substep(dt, step) parcels%vorticity(:, n) = parcels%vorticity(:, n) & + cb * dt * parcels%delta_vor(:, n) - call evolve_ellipsoid(parcels%B(:, n), parcels%strain(:, n), cb * dt) + call evolve_ellipsoid(parcels%B(:, n), parcels%int_strain(:, n), cb * dt) detB = parcels%B(1, n) * (parcels%B(4, n) * parcels%B(6, n) - parcels%B(5, n) ** 2) & - parcels%B(2, n) * (parcels%B(2, n) * parcels%B(6, n) - parcels%B(3, n) * parcels%B(5, n)) & @@ -191,7 +191,7 @@ subroutine ls_rk_substep(dt, step) do n = 1, n_parcels parcels%delta_pos(:, n) = ca * parcels%delta_pos(:, n) parcels%delta_vor(:, n) = ca * parcels%delta_vor(:, n) - parcels%strain(:, n) = ca * parcels%strain(:, n) + parcels%int_strain(:, n) = ca * parcels%int_strain(:, n) enddo !$omp end parallel do From 523edd63e99aaaee4dc1a1c74ebbd41b9dbe6ad5 Mon Sep 17 00:00:00 2001 From: sjboeing Date: Tue, 6 Feb 2024 10:27:12 +0000 Subject: [PATCH 47/62] Revert "Using one more iteration in the exponentiation" This reverts commit f178974a8b1a1e4d659d84668ba58dd5a3ff0282. --- src/3d/stepper/rk_utils.f90 | 7 +++---- 1 file changed, 3 insertions(+), 4 deletions(-) diff --git a/src/3d/stepper/rk_utils.f90 b/src/3d/stepper/rk_utils.f90 index d1fd63506..548acf418 100644 --- a/src/3d/stepper/rk_utils.f90 +++ b/src/3d/stepper/rk_utils.f90 @@ -50,11 +50,10 @@ subroutine evolve_ellipsoid(B, S, dt_sub) Smat(3, 3) = -(S(I_DUDX) + S(I_DVDY)) ! S33 Qmat = Imat + & - 0.125 * dt_sub * Smat + & - 0.0078125 * dt_sub * dt_sub * matmul(Smat, Smat) + & - (one / 3072.0) * dt_sub * dt_sub * dt_sub * matmul(Smat, matmul(Smat, Smat)) + 0.25 * dt_sub * Smat + & + 0.03125 * dt_sub * dt_sub * matmul(Smat, Smat) + & + (1.0 / 384.0) * dt_sub * dt_sub * dt_sub * matmul(Smat, matmul(Smat, Smat)) - Qmat = matmul(Qmat, Qmat) Qmat = matmul(Qmat, Qmat) Qmat = matmul(Qmat, Qmat) Bmat = matmul(Qmat, matmul(Bmat, transpose(Qmat))) From dcae3a923e90f8943a304abd97a0cd714849c30a Mon Sep 17 00:00:00 2001 From: sjboeing Date: Wed, 14 Feb 2024 11:55:02 +0000 Subject: [PATCH 48/62] Work on 2d version --- examples/straka.config | 2 +- examples/straka2d.nml | 2 +- src/2d/parcels/parcel_interpl.f90 | 4 +-- src/2d/stepper/ls_rk4.f90 | 29 ++++++------------- src/2d/stepper/rk4_utils.f90 | 47 ++++++++++++++++++++----------- 5 files changed, 43 insertions(+), 41 deletions(-) diff --git a/examples/straka.config b/examples/straka.config index 6ec443c40..5dde64a3b 100644 --- a/examples/straka.config +++ b/examples/straka.config @@ -9,7 +9,7 @@ output%field_freq = 20 - output%parcel_freq = 20 + output%parcel_freq = 100 output%parcel_stats_freq = 1 diff --git a/examples/straka2d.nml b/examples/straka2d.nml index 31a8fcdea..635c9a3a9 100644 --- a/examples/straka2d.nml +++ b/examples/straka2d.nml @@ -3,7 +3,7 @@ model = 'Straka' ncfname = 'straka2d.nc' - box%ncells = 128, 16 + box%ncells = 2048, 256 box%extent = 51200, 6400 box%origin = -25600.0, 0.0 diff --git a/src/2d/parcels/parcel_interpl.f90 b/src/2d/parcels/parcel_interpl.f90 index ac171e276..1f5d7783a 100644 --- a/src/2d/parcels/parcel_interpl.f90 +++ b/src/2d/parcels/parcel_interpl.f90 @@ -312,6 +312,7 @@ subroutine grid2par(vel, vor, vgrad, add) do n = 1, n_parcels vel(:, n) = zero vor(n) = zero + vgrad(:, n) = zero enddo !$omp end do !$omp end parallel @@ -322,6 +323,7 @@ subroutine grid2par(vel, vor, vgrad, add) do n = 1, n_parcels vel(:, n) = zero vor(n) = zero + vgrad(:, n) = zero enddo !$omp end do !$omp end parallel @@ -331,8 +333,6 @@ subroutine grid2par(vel, vor, vgrad, add) !$omp do private(n, p, l, points, weight, is, js, weights) do n = 1, n_parcels - vgrad(:, n) = zero - points = get_ellipse_points(parcels%position(:, n), & parcels%volume(n), & parcels%B(:, n)) diff --git a/src/2d/stepper/ls_rk4.f90 b/src/2d/stepper/ls_rk4.f90 index 66a50855d..bd072f1e9 100644 --- a/src/2d/stepper/ls_rk4.f90 +++ b/src/2d/stepper/ls_rk4.f90 @@ -6,7 +6,7 @@ module ls_rk4 use options, only : parcel use parcel_container use parcel_bc - use rk4_utils, only: get_B, get_time_step + use rk4_utils, only: evolve_ellipsoid, get_time_step use utils, only : write_step use parcel_interpl, only : par2grid, grid2par, grid2par_add use fields, only : velgradg, velog, vortg, vtend, tbuoyg @@ -23,7 +23,7 @@ module ls_rk4 double precision, allocatable, dimension(:, :) :: & delta_pos, & - strain, & ! strain at parcel location + int_strain, & ! strain at parcel location delta_b ! B matrix integration double precision, allocatable, dimension(:) :: & @@ -53,7 +53,7 @@ subroutine ls_rk4_alloc(num) allocate(delta_pos(2, num)) allocate(delta_vor(num)) - allocate(strain(4, num)) + allocate(int_strain(4, num)) allocate(delta_b(2, num)) end subroutine ls_rk4_alloc @@ -63,7 +63,7 @@ subroutine ls_rk4_dealloc deallocate(delta_pos) deallocate(delta_vor) - deallocate(strain) + deallocate(int_strain) deallocate(delta_b) end subroutine ls_rk4_dealloc @@ -90,7 +90,7 @@ subroutine ls_rk4_step(t) ! update the time step dt = get_time_step(t) - call grid2par(delta_pos, delta_vor, strain) + call grid2par(delta_pos, delta_vor, int_strain) call calculate_parcel_diagnostics(delta_pos) @@ -131,29 +131,16 @@ subroutine ls_rk4_substep(dt, step) if (step == 1) then call start_timer(rk4_timer) - !$omp parallel do default(shared) private(n) - do n = 1, n_parcels - delta_b(:, n) = get_B(parcels%B(:, n), strain(:, n), parcels%volume(n)) - enddo - !$omp end parallel do - call stop_timer(rk4_timer) else call vor2vel(vortg, velog, velgradg) call vorticity_tendency(tbuoyg, vtend) - call grid2par_add(delta_pos, delta_vor, strain) + call grid2par_add(delta_pos, delta_vor, int_strain) call start_timer(rk4_timer) - !$omp parallel do default(shared) private(n) - do n = 1, n_parcels - delta_b(:, n) = delta_b(:, n) & - + get_B(parcels%B(:, n), strain(:, n), parcels%volume(n)) - enddo - !$omp end parallel do - call stop_timer(rk4_timer) endif @@ -165,7 +152,7 @@ subroutine ls_rk4_substep(dt, step) + cb * dt * delta_pos(:, n) parcels%vorticity(n) = parcels%vorticity(n) + cb * dt * delta_vor(n) - parcels%B(:, n) = parcels%B(:, n) + cb * dt * delta_b(:, n) + call evolve_ellipsoid(parcels%B(:, n), int_strain(:, n), parcels%volume(n), cb * dt) enddo !$omp end parallel do @@ -181,7 +168,7 @@ subroutine ls_rk4_substep(dt, step) do n = 1, n_parcels delta_pos(:, n) = ca * delta_pos(:, n) delta_vor(n) = ca * delta_vor(n) - delta_b(:, n) = ca * delta_b(:, n) + int_strain(:, n) = ca * int_strain(:, n) enddo !$omp end parallel do diff --git a/src/2d/stepper/rk4_utils.f90 b/src/2d/stepper/rk4_utils.f90 index be3df451d..d6c1a37d9 100644 --- a/src/2d/stepper/rk4_utils.f90 +++ b/src/2d/stepper/rk4_utils.f90 @@ -8,27 +8,42 @@ module rk4_utils #endif implicit none + double precision, parameter :: Imat(2,2)=reshape((/one,zero,zero,one/), (/2,2/)) contains - ! Advance the B matrix. - ! @param[in] Bin are the B matrix components of the parcel - ! @param[in] S is the local velocity strain - ! @param[in] volume is the parcel volume - ! @returns the updated B matrix components (B11 and B12) in Bout - function get_B(Bin, S, volume) result(Bout) - double precision, intent(in) :: Bin(2) + subroutine evolve_ellipsoid(B, S, volume, dt_sub) + double precision, intent(inout) :: B(2) double precision, intent(in) :: S(4) double precision, intent(in) :: volume - double precision :: Bout(2) - - ! B11 = 2 * (dudx * B11 + dudy * B12) - Bout(1) = two * (S(1) * Bin(1) + S(2) * Bin(2)) - - ! B12 = dvdx * B11 + dudy * B22 - Bout(2) = S(3) * Bin(1) + S(2) * get_B22(Bin(1), Bin(2), volume) - - end function get_B + double precision, intent(in) :: dt_sub + double precision :: Bmat(2,2) + double precision :: Smat(2,2) + double precision :: Qmat(2,2) + + Bmat(1, 1) = B(1) ! B11 + Bmat(1, 2) = B(2) ! B12 + Bmat(2, 1) = B(1) ! B21 + Bmat(2, 2) = get_B22(B(1), B(2), volume) + + Smat(1, 1) = S(1) ! S11 + Smat(1, 2) = S(2) ! S12 + Smat(2, 1) = S(3) ! S21 + Smat(2, 2) = S(4) ! S33 + + Qmat = Imat + & + 0.25 * dt_sub * Smat + & + 0.03125 * dt_sub * dt_sub * matmul(Smat, Smat) + & + (1.0 / 384.0) * dt_sub * dt_sub * dt_sub * matmul(Smat, matmul(Smat, Smat)) + + Qmat = matmul(Qmat, Qmat) + Qmat = matmul(Qmat, Qmat) + Bmat = matmul(Qmat, matmul(Bmat, transpose(Qmat))) + + B(1) = Bmat(1, 1) + B(2) = Bmat(1, 2) + + end subroutine evolve_ellipsoid ! Estimate a suitable time step based on the velocity strain ! and buoyancy gradient. From 921363e9ab711d611e847247548d8a15f3f10376 Mon Sep 17 00:00:00 2001 From: sjboeing Date: Fri, 16 Feb 2024 08:15:44 +0000 Subject: [PATCH 49/62] Optimise 3D matrix exponentiation --- src/3d/parcels/parcel_interpl.f90 | 2 +- src/3d/stepper/rk_utils.f90 | 34 ++++++++++++++++++++++++------- 2 files changed, 28 insertions(+), 8 deletions(-) diff --git a/src/3d/parcels/parcel_interpl.f90 b/src/3d/parcels/parcel_interpl.f90 index 66814bb3e..576e9b4fc 100644 --- a/src/3d/parcels/parcel_interpl.f90 +++ b/src/3d/parcels/parcel_interpl.f90 @@ -461,7 +461,7 @@ subroutine grid2par(add) enddo do l = 1,8 parcels%int_strain(l, n) = parcels%int_strain(l, n) & - + point_weight_g2p * sum(weights * velgradg(ks:ks+1, js:js+1, is:is+1, l)) + + point_weight_g2p * sum(weights * velgradg(ks:ks+1, js:js+1, is:is+1, l)) enddo do l = 1,3 parcels%delta_vor(l, n) = parcels%delta_vor(l, n) & diff --git a/src/3d/stepper/rk_utils.f90 b/src/3d/stepper/rk_utils.f90 index 4752e73c8..df17438b8 100644 --- a/src/3d/stepper/rk_utils.f90 +++ b/src/3d/stepper/rk_utils.f90 @@ -3,7 +3,7 @@ module rk_utils use parcel_ellipsoid, only : get_B33, I_B11, I_B12, I_B13, I_B22, I_B23, I_B33 use fields, only : velgradg, tbuoyg, vortg, I_DUDX, I_DUDY, I_DUDZ, I_DVDX, I_DVDY, I_DVDZ, I_DWDX, I_DWDY, strain_mag use field_mpi, only : field_halo_fill_scalar - use constants, only : zero, one, two, f12 + use constants, only : zero, one, two, f12, f23 use parameters, only : nx, ny, nz, dxi, vcell use scherzinger, only : scherzinger_eigenvalues use mpi_layout, only : box @@ -17,7 +17,16 @@ module rk_utils #endif implicit none - double precision, parameter :: Imat(3,3)=reshape((/one,zero,zero,zero,one,zero,zero,zero,one/), (/3,3/)) + double precision, parameter :: Imat(3,3) = reshape((/one,zero,zero,zero,one,zero,zero,zero,one/), (/3,3/)) + double precision, parameter :: f177 = sqrt(177.0) + double precision, parameter :: xx1 = (one / 88.0) * (one + f177) * f23 + double precision, parameter :: xx2 = (one / 352.0) * (one + f177) * f23 + double precision, parameter :: xx4 = (-271.0 + 29.0 * f177) / (315.0 * f23) + double precision, parameter :: xx5 = (11.0 * (-1.0 + f177)) / (1260.0 * f23) + double precision, parameter :: xx6 = (11.0 * (-9.0 + f177)) / (5040.0 * f23) + double precision, parameter :: xx7 = (89.0 - f177) / (5040.0 * f23 * f23) + double precision, parameter :: yy2 = (one / 630.0) * (857.0 - 58.0 * f177) + contains @@ -27,6 +36,10 @@ subroutine evolve_ellipsoid(B, S, dt_sub) double precision, intent(in) :: dt_sub double precision :: Bmat(3,3) double precision :: Smat(3,3) + double precision :: Rmat(3,3) + double precision :: Rmat2(3,3) + double precision :: Rmat4(3,3) + double precision :: Rmat8(3,3) double precision :: Qmat(3,3) Bmat(1, 1) = B(I_B11) ! B11 @@ -49,11 +62,18 @@ subroutine evolve_ellipsoid(B, S, dt_sub) Smat(3, 2) = S(I_DWDY) ! S32 Smat(3, 3) = -(S(I_DUDX) + S(I_DVDY)) ! S33 - Qmat = Imat + & - 0.25 * dt_sub * Smat + & - 0.03125 * dt_sub * dt_sub * matmul(Smat, Smat) + & - (1.0 / 384.0) * dt_sub * dt_sub * dt_sub * matmul(Smat, matmul(Smat, Smat)) - + ! Bader, P., Blanes, S., & Casas, F. (2019). + ! Computing the matrix exponential with an optimized Taylor polynomial approximation. + ! Mathematics, 7(12), 1174. + ! Using 8th order Taylor with 2 ward steps + ! Possibly this is overkill and ward steps can be removed + ! This does not save much computation though + + Rmat = 0.25 * dt_sub * Smat + Rmat2 = matmul(Rmat, Rmat) + Rmat4 = matmul(Rmat2, xx1 * Rmat + xx2 * Rmat2) + Rmat8 = matmul(f23 * Rmat2 + Rmat4, xx4 * Imat + xx5 * Rmat + xx6 * Rmat2 + xx7 * Rmat4) + Qmat = Imat + Rmat + yy2 * Rmat2 + Rmat8 Qmat = matmul(Qmat, Qmat) Qmat = matmul(Qmat, Qmat) Bmat = matmul(Qmat, matmul(Bmat, transpose(Qmat))) From 50c498adb17facb5aba361aa9a663e8bee65709e Mon Sep 17 00:00:00 2001 From: sjboeing Date: Fri, 16 Feb 2024 09:19:26 +0000 Subject: [PATCH 50/62] Fancy matrix exp in 2d and 3d --- src/2d/stepper/rk4_utils.f90 | 23 +++++++++++++++++------ src/3d/stepper/rk_utils.f90 | 11 +---------- src/utils/constants.f90 | 9 +++++++++ 3 files changed, 27 insertions(+), 16 deletions(-) diff --git a/src/2d/stepper/rk4_utils.f90 b/src/2d/stepper/rk4_utils.f90 index d6c1a37d9..a26fc4112 100644 --- a/src/2d/stepper/rk4_utils.f90 +++ b/src/2d/stepper/rk4_utils.f90 @@ -1,7 +1,7 @@ module rk4_utils use parcel_ellipse, only : get_B22 use fields, only : velgradg, tbuoyg, vtend - use constants, only : zero, one, two, f12 + use constants, only : zero, one, two, f12, f23, xx1, xx2, xx4, xx5, xx6, xx7, yy2 use parameters, only : nx, nz, dxi #ifdef ENABLE_VERBOSE use options, only : output @@ -20,6 +20,10 @@ subroutine evolve_ellipsoid(B, S, volume, dt_sub) double precision :: Bmat(2,2) double precision :: Smat(2,2) double precision :: Qmat(2,2) + double precision :: Rmat(2,2) + double precision :: Rmat2(2,2) + double precision :: Rmat4(2,2) + double precision :: Rmat8(2,2) Bmat(1, 1) = B(1) ! B11 Bmat(1, 2) = B(2) ! B12 @@ -31,11 +35,18 @@ subroutine evolve_ellipsoid(B, S, volume, dt_sub) Smat(2, 1) = S(3) ! S21 Smat(2, 2) = S(4) ! S33 - Qmat = Imat + & - 0.25 * dt_sub * Smat + & - 0.03125 * dt_sub * dt_sub * matmul(Smat, Smat) + & - (1.0 / 384.0) * dt_sub * dt_sub * dt_sub * matmul(Smat, matmul(Smat, Smat)) - + ! Bader, P., Blanes, S., & Casas, F. (2019). + ! Computing the matrix exponential with an optimized Taylor polynomial approximation. + ! Mathematics, 7(12), 1174. + ! Using 8th order Taylor with 2 ward steps + ! Possibly this is overkill and ward steps can be removed + ! This does not save much computation though + + Rmat = 0.25 * dt_sub * Smat + Rmat2 = matmul(Rmat, Rmat) + Rmat4 = matmul(Rmat2, xx1 * Rmat + xx2 * Rmat2) + Rmat8 = matmul(f23 * Rmat2 + Rmat4, xx4 * Imat + xx5 * Rmat + xx6 * Rmat2 + xx7 * Rmat4) + Qmat = Imat + Rmat + yy2 * Rmat2 + Rmat8 Qmat = matmul(Qmat, Qmat) Qmat = matmul(Qmat, Qmat) Bmat = matmul(Qmat, matmul(Bmat, transpose(Qmat))) diff --git a/src/3d/stepper/rk_utils.f90 b/src/3d/stepper/rk_utils.f90 index df17438b8..b7fd2e6ea 100644 --- a/src/3d/stepper/rk_utils.f90 +++ b/src/3d/stepper/rk_utils.f90 @@ -3,7 +3,7 @@ module rk_utils use parcel_ellipsoid, only : get_B33, I_B11, I_B12, I_B13, I_B22, I_B23, I_B33 use fields, only : velgradg, tbuoyg, vortg, I_DUDX, I_DUDY, I_DUDZ, I_DVDX, I_DVDY, I_DVDZ, I_DWDX, I_DWDY, strain_mag use field_mpi, only : field_halo_fill_scalar - use constants, only : zero, one, two, f12, f23 + use constants, only : zero, one, two, f12, f23, xx1, xx2, xx4, xx5, xx6, xx7, yy2 use parameters, only : nx, ny, nz, dxi, vcell use scherzinger, only : scherzinger_eigenvalues use mpi_layout, only : box @@ -18,15 +18,6 @@ module rk_utils implicit none double precision, parameter :: Imat(3,3) = reshape((/one,zero,zero,zero,one,zero,zero,zero,one/), (/3,3/)) - double precision, parameter :: f177 = sqrt(177.0) - double precision, parameter :: xx1 = (one / 88.0) * (one + f177) * f23 - double precision, parameter :: xx2 = (one / 352.0) * (one + f177) * f23 - double precision, parameter :: xx4 = (-271.0 + 29.0 * f177) / (315.0 * f23) - double precision, parameter :: xx5 = (11.0 * (-1.0 + f177)) / (1260.0 * f23) - double precision, parameter :: xx6 = (11.0 * (-9.0 + f177)) / (5040.0 * f23) - double precision, parameter :: xx7 = (89.0 - f177) / (5040.0 * f23 * f23) - double precision, parameter :: yy2 = (one / 630.0) * (857.0 - 58.0 * f177) - contains diff --git a/src/utils/constants.f90 b/src/utils/constants.f90 index 5cfd74699..0db387ffe 100644 --- a/src/utils/constants.f90 +++ b/src/utils/constants.f90 @@ -45,4 +45,13 @@ module constants double precision, parameter :: rad2deg = 180.0d0 * fpi double precision, parameter :: deg2rad = one / rad2deg + double precision, parameter :: sq177 = sqrt(177.0) + double precision, parameter :: xx1 = (one / 88.0) * (one + sq177) * f23 + double precision, parameter :: xx2 = (one / 352.0) * (one + sq177) * f23 + double precision, parameter :: xx4 = (-271.0 + 29.0 * sq177) / (315.0 * f23) + double precision, parameter :: xx5 = (11.0 * (-1.0 + sq177)) / (1260.0 * f23) + double precision, parameter :: xx6 = (11.0 * (-9.0 + sq177)) / (5040.0 * f23) + double precision, parameter :: xx7 = (89.0 - sq177) / (5040.0 * f23 * f23) + double precision, parameter :: yy2 = (one / 630.0) * (857.0 - 58.0 * sq177) + end module From a727f96997b34dc71b8c8b854f7e0cfaa13aaee4 Mon Sep 17 00:00:00 2001 From: sjboeing Date: Wed, 21 Feb 2024 17:11:01 +0000 Subject: [PATCH 51/62] Correct B only at end of time step --- src/3d/stepper/ls_rk.f90 | 33 +++++++++++++++++---------------- 1 file changed, 17 insertions(+), 16 deletions(-) diff --git a/src/3d/stepper/ls_rk.f90 b/src/3d/stepper/ls_rk.f90 index 71a611da6..d94c2505c 100644 --- a/src/3d/stepper/ls_rk.f90 +++ b/src/3d/stepper/ls_rk.f90 @@ -178,13 +178,12 @@ subroutine ls_rk_substep(dt, step) enddo !$omp end parallel do - call stop_timer(rk_timer) endif call start_timer(rk_timer) - !$omp parallel do default(shared) private(n,detB,factor) + !$omp parallel do default(shared) private(n) do n = 1, n_parcels parcels%position(:, n) = parcels%position(:, n) & + cb * dt * parcels%delta_pos(:, n) @@ -194,28 +193,30 @@ subroutine ls_rk_substep(dt, step) parcels%B(:, n) = parcels%B(:, n) & + cb * dt * parcels%delta_b(:, n) - - ! normalize B matrix - detB = parcels%B(1, n) * (parcels%B(4, n) * parcels%B(6, n) - parcels%B(5, n) ** 2) & - - parcels%B(2, n) * (parcels%B(2, n) * parcels%B(6, n) - parcels%B(3, n) * parcels%B(5, n)) & - + parcels%B(3, n) * (parcels%B(2, n) * parcels%B(5, n) - parcels%B(3, n) * parcels%B(4, n)) - - factor = (get_abc(parcels%volume(n)) ** 2 / detB) ** f13 - - parcels%delta_b(:, n) = parcels%delta_b(:, n) + parcels%B(:, n) * (factor - one)/(cb * dt) - parcels%B(:, n) = parcels%B(:, n) * factor enddo !$omp end parallel do - call stop_timer(rk_timer) - if (step == n_stages) then + !$omp parallel do default(shared) private(n,detB,factor) + do n = 1, n_parcels + ! normalize B matrix in final substep + detB = parcels%B(1, n) * (parcels%B(4, n) * parcels%B(6, n) - parcels%B(5, n) ** 2) & + - parcels%B(2, n) * (parcels%B(2, n) * parcels%B(6, n) - parcels%B(3, n) * parcels%B(5, n)) & + + parcels%B(3, n) * (parcels%B(2, n) * parcels%B(5, n) - parcels%B(3, n) * parcels%B(4, n)) + + factor = (get_abc(parcels%volume(n)) ** 2 / detB) ** f13 + + parcels%B(:, n) = parcels%B(:, n) * factor + enddo + !$omp end parallel do + call parcel_communicate + + call stop_timer(rk_timer) + return endif - call start_timer(rk_timer) - !$omp parallel do default(shared) private(n) do n = 1, n_parcels parcels%delta_pos(:, n) = ca * parcels%delta_pos(:, n) From 231ebba4fd1507d9f4a0d2e98871f6e2d980afa7 Mon Sep 17 00:00:00 2001 From: sjboeing Date: Thu, 22 Feb 2024 10:17:07 +0000 Subject: [PATCH 52/62] Fix parcel split --- unit-tests/3d/test_mpi_parcel_split.f90 | 11 ++++++----- 1 file changed, 6 insertions(+), 5 deletions(-) diff --git a/unit-tests/3d/test_mpi_parcel_split.f90 b/unit-tests/3d/test_mpi_parcel_split.f90 index c75c648d3..248e2be90 100644 --- a/unit-tests/3d/test_mpi_parcel_split.f90 +++ b/unit-tests/3d/test_mpi_parcel_split.f90 @@ -7,7 +7,7 @@ ! domain of *this* process. ! ============================================================================= program test_mpi_parcel_split - use constants, only : zero, one, f18, f14, f12, f34, fpi, pi, four + use constants, only : zero, one, four, f14, f12, f23, f34, fpi, pi use unit_test use mpi_environment use mpi_layout @@ -57,12 +57,13 @@ program test_mpi_parcel_split ! volume per parcel is f12 * vcell ! f12 * vcell = four / three * abc * pi --> abc = f34 * f12 * vcell * fpi abc = f34 * f12 * vcell * fpi - a2 = f34 * abc - b2 = f18 * abc + a2 = 16.0d0 * abc**f23 + b2 = f14 * abc**f23 c2 = b2 - abc = dsqrt(a2) - call set_amax(abc) + ! Set amax larger than a2 + ! Force split by aspect ratio set above, rather than amax + call set_amax(1.1d0*sqrt(a2)) ! place parcels in the last interior cells in the west i = box%lo(1) From cf0b47af501ba6ec8e418176ae2cc6e7736e71e4 Mon Sep 17 00:00:00 2001 From: sjboeing Date: Sun, 18 Feb 2024 19:42:02 +0000 Subject: [PATCH 53/62] Temporary fixes for debug compilation --- src/3d/parcels/parcel_damping.f90 | 7 ++++++- src/3d/parcels/parcel_interpl.f90 | 8 ++++++++ 2 files changed, 14 insertions(+), 1 deletion(-) diff --git a/src/3d/parcels/parcel_damping.f90 b/src/3d/parcels/parcel_damping.f90 index c442ed98a..4589c7106 100644 --- a/src/3d/parcels/parcel_damping.f90 +++ b/src/3d/parcels/parcel_damping.f90 @@ -70,7 +70,7 @@ end subroutine parcel_damp ! @pre subroutine perturbation_damping(dt, l_reuse) double precision, intent(in) :: dt - logical, intent(in) :: l_reuse + logical :: l_reuse integer :: n, p, l double precision :: points(3, n_points_p2g) double precision :: pvol @@ -84,6 +84,11 @@ subroutine perturbation_damping(dt, l_reuse) call start_timer(damping_timer) + ! This is only here to allow debug compilation +#ifdef ENABLE_P2G_1POINT + l_reuse=l_reuse +#endif + !$omp parallel default(shared) !$omp do private(n, p, l, points, pvol, weight) & #ifndef ENABLE_DRY_MODE diff --git a/src/3d/parcels/parcel_interpl.f90 b/src/3d/parcels/parcel_interpl.f90 index e4dd83768..b2131f5d3 100644 --- a/src/3d/parcels/parcel_interpl.f90 +++ b/src/3d/parcels/parcel_interpl.f90 @@ -148,6 +148,7 @@ end subroutine vol2grid ! - nparg, that is the number of parcels per grid cell ! - nsparg, that is the number of small parcels per grid cell ! @pre The parcel must be assigned to the correct MPI process. + subroutine par2grid(l_reuse) logical, optional :: l_reuse double precision :: points(3, n_points_p2g) @@ -159,6 +160,13 @@ subroutine par2grid(l_reuse) call start_timer(par2grid_timer) + ! This is only here to allow debug compilation +#ifdef ENABLE_P2G_1POINT + if(present(l_reuse)) then + l_reuse=l_reuse + endif +#endif + vortg = zero volg = zero nparg = zero From b07df4833e0680ea601c18e05db11a32a193f8f9 Mon Sep 17 00:00:00 2001 From: sjboeing Date: Thu, 22 Feb 2024 12:20:23 +0000 Subject: [PATCH 54/62] Fix debug compilation in 1-point p2g and g2p mode --- src/3d/parcels/parcel_damping.f90 | 5 +++-- src/3d/parcels/parcel_interpl.f90 | 5 +++-- 2 files changed, 6 insertions(+), 4 deletions(-) diff --git a/src/3d/parcels/parcel_damping.f90 b/src/3d/parcels/parcel_damping.f90 index 4589c7106..33da4dda0 100644 --- a/src/3d/parcels/parcel_damping.f90 +++ b/src/3d/parcels/parcel_damping.f90 @@ -85,8 +85,9 @@ subroutine perturbation_damping(dt, l_reuse) call start_timer(damping_timer) ! This is only here to allow debug compilation -#ifdef ENABLE_P2G_1POINT - l_reuse=l_reuse + ! with a warning for unused variables +#if defined (ENABLE_P2G_1POINT) && !defined (NDEBUG) + l_reuse=.false. #endif !$omp parallel default(shared) diff --git a/src/3d/parcels/parcel_interpl.f90 b/src/3d/parcels/parcel_interpl.f90 index b2131f5d3..7dce1ce20 100644 --- a/src/3d/parcels/parcel_interpl.f90 +++ b/src/3d/parcels/parcel_interpl.f90 @@ -161,9 +161,10 @@ subroutine par2grid(l_reuse) call start_timer(par2grid_timer) ! This is only here to allow debug compilation -#ifdef ENABLE_P2G_1POINT + ! with a warning for unused variables +#if defined (ENABLE_P2G_1POINT) && !defined (NDEBUG) if(present(l_reuse)) then - l_reuse=l_reuse + l_reuse=.false. endif #endif From 16b0567cc190ca584d998de870ae1b0c538b3773 Mon Sep 17 00:00:00 2001 From: sjboeing Date: Thu, 22 Feb 2024 17:45:33 +0000 Subject: [PATCH 55/62] Fix unit test issue with renaming parcel strain --- unit-tests/3d/test_mpi_grid2par.f90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/unit-tests/3d/test_mpi_grid2par.f90 b/unit-tests/3d/test_mpi_grid2par.f90 index 4b2d7f873..2b1572a20 100644 --- a/unit-tests/3d/test_mpi_grid2par.f90 +++ b/unit-tests/3d/test_mpi_grid2par.f90 @@ -103,7 +103,7 @@ program test_mpi_grid2par enddo do l = 1, 5 - error = max(error, maxval(dabs(parcels%strain(l, 1:n_parcels) - dble(l)))) + error = max(error, maxval(dabs(parcels%int_strain(l, 1:n_parcels) - dble(l)))) enddo call mpi_blocking_reduce(error, MPI_MAX, world) From 859afac6b032435706e9c6316070d2bb42f9feff Mon Sep 17 00:00:00 2001 From: sjboeing Date: Mon, 5 Aug 2024 12:43:04 +0100 Subject: [PATCH 56/62] Remove unneeded calls --- src/3d/stepper/ls_rk.f90 | 23 +++++++---------------- 1 file changed, 7 insertions(+), 16 deletions(-) diff --git a/src/3d/stepper/ls_rk.f90 b/src/3d/stepper/ls_rk.f90 index 958910883..fc8dd7df6 100644 --- a/src/3d/stepper/ls_rk.f90 +++ b/src/3d/stepper/ls_rk.f90 @@ -26,7 +26,6 @@ module ls_rk private integer, parameter :: dp=kind(zero) ! double precision - integer :: rk_timer ! fourth order RK coefficients: @@ -131,9 +130,9 @@ subroutine ls_rk_step(t) call parcel_damp(dt) - ! we need to subtract 14 calls since we start and stop - ! the timer multiple times which increments n_calls - timings(rk_timer)%n_calls = timings(rk_timer)%n_calls - (3 * n_stages - 1) + ! we need to subtract 1 call since we start and stop + ! the timer again when we apply the reflective bc + timings(rk_timer)%n_calls = timings(rk_timer)%n_calls - 1 t = t + dt end subroutine ls_rk_step @@ -152,20 +151,12 @@ subroutine ls_rk_substep(dt, step) ca = captr(step) cb = cbptr(step) - if (step == 1) then - call start_timer(rk_timer) - - call stop_timer(rk_timer) - else + if (step > 1) then call vor2vel call vorticity_tendency call grid2par(add=.true.) - - call start_timer(rk_timer) - - call stop_timer(rk_timer) endif call start_timer(rk_timer) @@ -177,14 +168,14 @@ subroutine ls_rk_substep(dt, step) parcels%vorticity(:, n) = parcels%vorticity(:, n) & + cb * dt * parcels%delta_vor(:, n) - + call evolve_ellipsoid(parcels%B(:, n), parcels%int_strain(:, n), cb * dt) - + enddo !$omp end parallel do if (step == n_stages) then - !$omp parallel do default(shared) private(n,detB,factor) + !$omp parallel do default(shared) private(n, detB, factor) do n = 1, n_parcels ! normalize B matrix in final substep detB = parcels%B(1, n) * (parcels%B(4, n) * parcels%B(6, n) - parcels%B(5, n) ** 2) & From a88ded1bfa5911e83ba5dd7e90ad00bf3d2cd3d8 Mon Sep 17 00:00:00 2001 From: sjboeing Date: Mon, 5 Aug 2024 12:43:43 +0100 Subject: [PATCH 57/62] Minor optimisations in rk utils --- src/3d/stepper/rk_utils.f90 | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/src/3d/stepper/rk_utils.f90 b/src/3d/stepper/rk_utils.f90 index 8f0cbf060..ba9c4b50e 100644 --- a/src/3d/stepper/rk_utils.f90 +++ b/src/3d/stepper/rk_utils.f90 @@ -3,7 +3,7 @@ module rk_utils use parcel_ellipsoid, only : get_B33, I_B11, I_B12, I_B13, I_B22, I_B23, I_B33 use fields, only : velgradg, tbuoyg, vortg, I_DUDX, I_DUDY, I_DUDZ, I_DVDX, I_DVDY, I_DVDZ, I_DWDX, I_DWDY, strain_mag use field_mpi, only : field_halo_fill_scalar - use constants, only : zero, one, two, f12, f23, xx1, xx2, xx4, xx5, xx6, xx7, yy2 + use constants, only : zero, one, two, f12, f14, f23, xx1, xx2, xx4, xx5, xx6, xx7, yy2 use parameters, only : nx, ny, nz, dxi, vcell use scherzinger, only : scherzinger_eigenvalues use mpi_layout, only : box @@ -36,11 +36,11 @@ subroutine evolve_ellipsoid(B, S, dt_sub) Bmat(1, 1) = B(I_B11) ! B11 Bmat(1, 2) = B(I_B12) ! B12 Bmat(1, 3) = B(I_B13) ! B13 - Bmat(2, 1) = B(I_B12) ! B21 + Bmat(2, 1) = Bmat(1, 2) ! B21 Bmat(2, 2) = B(I_B22) ! B22 Bmat(2, 3) = B(I_B23) ! B23 - Bmat(3, 1) = B(I_B13) ! B31 - Bmat(3, 2) = B(I_B23) ! B32 + Bmat(3, 1) = Bmat(1, 3) ! B31 + Bmat(3, 2) = Bmat(2, 3) ! B32 Bmat(3, 3) = B(I_B33) ! B33 Smat(1, 1) = S(I_DUDX) ! S11 @@ -51,7 +51,7 @@ subroutine evolve_ellipsoid(B, S, dt_sub) Smat(2, 3) = S(I_DVDZ) ! S23 Smat(3, 1) = S(I_DWDX) ! S31 Smat(3, 2) = S(I_DWDY) ! S32 - Smat(3, 3) = -(S(I_DUDX) + S(I_DVDY)) ! S33 + Smat(3, 3) = -(Smat(1, 1) + Smat(2, 2)) ! S33 ! Bader, P., Blanes, S., & Casas, F. (2019). ! Computing the matrix exponential with an optimized Taylor polynomial approximation. @@ -60,7 +60,7 @@ subroutine evolve_ellipsoid(B, S, dt_sub) ! Possibly this is overkill and ward steps can be removed ! This does not save much computation though - Rmat = 0.25 * dt_sub * Smat + Rmat = f14 * dt_sub * Smat Rmat2 = matmul(Rmat, Rmat) Rmat4 = matmul(Rmat2, xx1 * Rmat + xx2 * Rmat2) Rmat8 = matmul(f23 * Rmat2 + Rmat4, xx4 * Imat + xx5 * Rmat + xx6 * Rmat2 + xx7 * Rmat4) From ee03e03e4ca594509b143e3e4fe6665e0e9c24ce Mon Sep 17 00:00:00 2001 From: sjboeing Date: Mon, 5 Aug 2024 12:44:10 +0100 Subject: [PATCH 58/62] Write B33 in netCDF --- src/3d/parcels/parcel_netcdf.f90 | 33 ++++++++++++++++++++++++-------- 1 file changed, 25 insertions(+), 8 deletions(-) diff --git a/src/3d/parcels/parcel_netcdf.f90 b/src/3d/parcels/parcel_netcdf.f90 index 442790c28..686b1776a 100644 --- a/src/3d/parcels/parcel_netcdf.f90 +++ b/src/3d/parcels/parcel_netcdf.f90 @@ -55,17 +55,18 @@ module parcel_netcdf , NC_B12 = 15 & , NC_B13 = 16 & , NC_B22 = 17 & - , NC_B23 = 18 + , NC_B23 = 18 & + , NC_B33 = 19 logical :: l_first_write = .true. logical :: l_unable = .false. #ifndef ENABLE_DRY_MODE - integer, parameter :: NC_HUM = 19 + integer, parameter :: NC_HUM = 20 #ifdef ENABLE_LABELS - integer, parameter :: NC_LABEL = 20 & - , NC_DILUTION = 21 + integer, parameter :: NC_LABEL = 21 & + , NC_DILUTION = 22 type(netcdf_info) :: nc_dset(NC_DILUTION) #else @@ -75,12 +76,12 @@ module parcel_netcdf #else #ifdef ENABLE_LABELS - integer, parameter :: NC_LABEL = 19 & - , NC_DILUTION = 20 + integer, parameter :: NC_LABEL = 20 & + , NC_DILUTION = 21 type(netcdf_info) :: nc_dset(NC_DILUTION) #else - type(netcdf_info) :: nc_dset(NC_B23) + type(netcdf_info) :: nc_dset(NC_B33) #endif #endif @@ -295,6 +296,7 @@ subroutine write_netcdf_parcels(t) call write_parcel_attribute(NC_B13, parcels%B(3, :), start, cnt) call write_parcel_attribute(NC_B22, parcels%B(4, :), start, cnt) call write_parcel_attribute(NC_B23, parcels%B(5, :), start, cnt) + call write_parcel_attribute(NC_B33, parcels%B(6, :), start, cnt) call write_parcel_attribute(NC_VOL, parcels%volume, start, cnt) @@ -603,6 +605,13 @@ subroutine read_chunk(first, last, pfirst) "The parcel shape component B23 must be present! Exiting.") endif + if (has_dataset(ncid, 'B33')) then + call read_netcdf_dataset(ncid, 'B33', parcels%B(6, pfirst:plast), start, cnt) + else + call mpi_exit_on_error(& + "The parcel shape component B33 must be present! Exiting.") + endif + if (has_dataset(ncid, 'x_position')) then call read_netcdf_dataset(ncid, 'x_position', & parcels%position(1, pfirst:plast), start, cnt) @@ -719,6 +728,7 @@ subroutine set_netcdf_parcel_output nc_dset(NC_B13)%l_enabled = .true. nc_dset(NC_B22)%l_enabled = .true. nc_dset(NC_B23)%l_enabled = .true. + nc_dset(NC_B33)%l_enabled = .true. else ! check individual fields do n = 1, size(nc_dset) @@ -765,6 +775,7 @@ subroutine set_netcdf_parcel_output nc_dset(NC_B13)%l_enabled = .true. nc_dset(NC_B22)%l_enabled = .true. nc_dset(NC_B23)%l_enabled = .true. + nc_dset(NC_B33)%l_enabled = .true. endif #ifdef ENABLE_VERBOSE @@ -787,7 +798,7 @@ subroutine set_netcdf_parcel_output l_enabled_restart = (l_enabled_restart .and. nc_dset(n)%l_enabled) enddo - do n = NC_B11, NC_B23 + do n = NC_B11, NC_B33 l_enabled_restart = (l_enabled_restart .and. nc_dset(n)%l_enabled) enddo l_enabled_restart = (l_enabled_restart .and. nc_dset(NC_VOL)%l_enabled) @@ -893,6 +904,12 @@ subroutine set_netcdf_parcel_info unit='m^2', & dtype=NF90_DOUBLE) + nc_dset(NC_B33) = netcdf_info(name='B33', & + long_name='B33 element of shape matrix', & + std_name='', & + unit='m^2', & + dtype=NF90_DOUBLE) + nc_dset(NC_VOL) = netcdf_info(name='volume', & long_name='parcel volume', & std_name='', & From 6cdbeb1ff20680d069b44807f152c6797b3ea28a Mon Sep 17 00:00:00 2001 From: sjboeing Date: Thu, 8 Aug 2024 00:14:52 +0100 Subject: [PATCH 59/62] Rename matrix exponential constants --- src/2d/stepper/rk4_utils.f90 | 18 +++++++++++++----- src/3d/stepper/rk_utils.f90 | 9 +++++---- src/utils/constants.f90 | 22 ++++++++++++++-------- 3 files changed, 32 insertions(+), 17 deletions(-) diff --git a/src/2d/stepper/rk4_utils.f90 b/src/2d/stepper/rk4_utils.f90 index a26fc4112..971d2f2e3 100644 --- a/src/2d/stepper/rk4_utils.f90 +++ b/src/2d/stepper/rk4_utils.f90 @@ -1,7 +1,8 @@ module rk4_utils use parcel_ellipse, only : get_B22 use fields, only : velgradg, tbuoyg, vtend - use constants, only : zero, one, two, f12, f23, xx1, xx2, xx4, xx5, xx6, xx7, yy2 + use constants, only : zero, one, two, f12, f14, f23, c_matexp_x1, c_matexp_x2, c_matexp_x4, & + c_matexp_x5, c_matexp_x6, c_matexp_x7, c_matexp_y2 use parameters, only : nx, nz, dxi #ifdef ENABLE_VERBOSE use options, only : output @@ -42,11 +43,18 @@ subroutine evolve_ellipsoid(B, S, volume, dt_sub) ! Possibly this is overkill and ward steps can be removed ! This does not save much computation though - Rmat = 0.25 * dt_sub * Smat + ! Bader, P., Blanes, S., & Casas, F. (2019). + ! Computing the matrix exponential with an optimized Taylor polynomial approximation. + ! Mathematics, 7(12), 1174. + ! Using 8th order Taylor with 2 ward steps + ! Possibly this is overkill and ward steps can be removed + ! This does not save much computation though + + Rmat = f14 * dt_sub * Smat Rmat2 = matmul(Rmat, Rmat) - Rmat4 = matmul(Rmat2, xx1 * Rmat + xx2 * Rmat2) - Rmat8 = matmul(f23 * Rmat2 + Rmat4, xx4 * Imat + xx5 * Rmat + xx6 * Rmat2 + xx7 * Rmat4) - Qmat = Imat + Rmat + yy2 * Rmat2 + Rmat8 + Rmat4 = matmul(Rmat2, c_matexp_x1 * Rmat + c_matexp_x2 * Rmat2) + Rmat8 = matmul(f23 * Rmat2 + Rmat4, c_matexp_x4 * Imat + c_matexp_x5 * Rmat + c_matexp_x6 * Rmat2 + c_matexp_x7 * Rmat4) + Qmat = Imat + Rmat + c_matexp_y2 * Rmat2 + Rmat8 Qmat = matmul(Qmat, Qmat) Qmat = matmul(Qmat, Qmat) Bmat = matmul(Qmat, matmul(Bmat, transpose(Qmat))) diff --git a/src/3d/stepper/rk_utils.f90 b/src/3d/stepper/rk_utils.f90 index 4485af877..efed42c3c 100644 --- a/src/3d/stepper/rk_utils.f90 +++ b/src/3d/stepper/rk_utils.f90 @@ -3,7 +3,8 @@ module rk_utils use parcel_ellipsoid, only : get_B33, I_B11, I_B12, I_B13, I_B22, I_B23, I_B33 use fields, only : velgradg, tbuoyg, vortg, I_DUDX, I_DUDY, I_DUDZ, I_DVDX, I_DVDY, I_DVDZ, I_DWDX, I_DWDY, strain_mag use field_mpi, only : field_halo_fill_scalar - use constants, only : zero, one, two, f12, f14, f23, xx1, xx2, xx4, xx5, xx6, xx7, yy2 + use constants, only : zero, one, two, f12, f14, f23, c_matexp_x1, c_matexp_x2, c_matexp_x4, & + c_matexp_x5, c_matexp_x6, c_matexp_x7, c_matexp_y2 use parameters, only : nx, ny, nz, dxi, vcell use scherzinger, only : scherzinger_eigenvalues use mpi_layout, only : box @@ -62,9 +63,9 @@ subroutine evolve_ellipsoid(B, S, dt_sub) Rmat = f14 * dt_sub * Smat Rmat2 = matmul(Rmat, Rmat) - Rmat4 = matmul(Rmat2, xx1 * Rmat + xx2 * Rmat2) - Rmat8 = matmul(f23 * Rmat2 + Rmat4, xx4 * Imat + xx5 * Rmat + xx6 * Rmat2 + xx7 * Rmat4) - Qmat = Imat + Rmat + yy2 * Rmat2 + Rmat8 + Rmat4 = matmul(Rmat2, c_matexp_x1 * Rmat + c_matexp_x2 * Rmat2) + Rmat8 = matmul(f23 * Rmat2 + Rmat4, c_matexp_x4 * Imat + c_matexp_x5 * Rmat + c_matexp_x6 * Rmat2 + c_matexp_x7 * Rmat4) + Qmat = Imat + Rmat + c_matexp_y2 * Rmat2 + Rmat8 Qmat = matmul(Qmat, Qmat) Qmat = matmul(Qmat, Qmat) Bmat = matmul(Qmat, matmul(Bmat, transpose(Qmat))) diff --git a/src/utils/constants.f90 b/src/utils/constants.f90 index 0db387ffe..e92256da2 100644 --- a/src/utils/constants.f90 +++ b/src/utils/constants.f90 @@ -45,13 +45,19 @@ module constants double precision, parameter :: rad2deg = 180.0d0 * fpi double precision, parameter :: deg2rad = one / rad2deg - double precision, parameter :: sq177 = sqrt(177.0) - double precision, parameter :: xx1 = (one / 88.0) * (one + sq177) * f23 - double precision, parameter :: xx2 = (one / 352.0) * (one + sq177) * f23 - double precision, parameter :: xx4 = (-271.0 + 29.0 * sq177) / (315.0 * f23) - double precision, parameter :: xx5 = (11.0 * (-1.0 + sq177)) / (1260.0 * f23) - double precision, parameter :: xx6 = (11.0 * (-9.0 + sq177)) / (5040.0 * f23) - double precision, parameter :: xx7 = (89.0 - sq177) / (5040.0 * f23 * f23) - double precision, parameter :: yy2 = (one / 630.0) * (857.0 - 58.0 * sq177) + ! Bader, P., Blanes, S., & Casas, F. (2019). + ! Computing the matrix exponential with an optimized Taylor polynomial approximation. + ! Mathematics, 7(12), 1174. + ! Coefficients needed for 8th order Taylor series + ! Naming of coefficients follows this article, with prefix 'c_matexp' + + double precision, parameter :: sq177 = sqrt(177.d0) + double precision, parameter :: c_matexp_x1 = (one / 88.d0) * (one + sq177) * f23 + double precision, parameter :: c_matexp_x2 = (one / 352.d0) * (one + sq177) * f23 + double precision, parameter :: c_matexp_x4 = (-271.d0 + 29.d0 * sq177) / (315.d0 * f23) + double precision, parameter :: c_matexp_x5 = (11.d0 * (-1.d0 + sq177)) / (1260.d0 * f23) + double precision, parameter :: c_matexp_x6 = (11.d0 * (-9.d0 + sq177)) / (5040.d0 * f23) + double precision, parameter :: c_matexp_x7 = (89.d0 - sq177) / (5040.d0 * f23 * f23) + double precision, parameter :: c_matexp_y2 = (one / 630.d0) * (857.d0 - 58.d0 * sq177) end module From 81489cde466cea2ca4faf103c19e4027f847d7a6 Mon Sep 17 00:00:00 2001 From: sjboeing Date: Thu, 26 Sep 2024 16:48:10 +0100 Subject: [PATCH 60/62] Fixes in debugging --- src/2d/fields/fields.f90 | 3 +++ src/3d/parcels/parcel_netcdf.f90 | 6 +++++- 2 files changed, 8 insertions(+), 1 deletion(-) diff --git a/src/2d/fields/fields.f90 b/src/2d/fields/fields.f90 index 06c86190c..61b417738 100644 --- a/src/2d/fields/fields.f90 +++ b/src/2d/fields/fields.f90 @@ -80,6 +80,9 @@ subroutine field_default #endif nparg = zero nsparg = zero +#ifndef NDEBUG + sym_volg = zero +#endif end subroutine ! Get the lower index of the cell the parcel is in. diff --git a/src/3d/parcels/parcel_netcdf.f90 b/src/3d/parcels/parcel_netcdf.f90 index 686b1776a..6a95bdad0 100644 --- a/src/3d/parcels/parcel_netcdf.f90 +++ b/src/3d/parcels/parcel_netcdf.f90 @@ -556,7 +556,11 @@ subroutine read_chunk(first, last, pfirst) integer, intent(in) :: first, last, pfirst logical :: l_valid = .false. integer :: cnt(2), start(2) - integer :: num, plast, n + integer :: num, plast +#ifdef ENABLE_LABELS + integer :: n +#endif + num = last - first + 1 plast = pfirst + num - 1 From 049ffe54c80b7bb94b9a6b533c3859e7c69c2c65 Mon Sep 17 00:00:00 2001 From: sjboeing Date: Thu, 26 Sep 2024 16:48:48 +0100 Subject: [PATCH 61/62] Comments on parcel container (3D) --- src/3d/parcels/parcel_container.f90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/3d/parcels/parcel_container.f90 b/src/3d/parcels/parcel_container.f90 index 188606068..8e54ea83d 100644 --- a/src/3d/parcels/parcel_container.f90 +++ b/src/3d/parcels/parcel_container.f90 @@ -77,8 +77,8 @@ module parcel_container position, & vorticity, & B ! B matrix entries; ordering: - ! B(:, 1) = B11, B(:, 2) = B12, B(:, 3) = B13 - ! B(:, 4) = B22, B(:, 5) = B23, B(:, 6) = B33 + ! B(1, :) = B11, B(2, :) = B12, B(3, :) = B13 + ! B(4, :) = B22, B(5, :) = B23, B(6, :) = B33 double precision, allocatable, dimension(:) :: & volume, & From 4ccf7f19abc59d06464c077ba14061fab57869d4 Mon Sep 17 00:00:00 2001 From: sjboeing Date: Thu, 26 Sep 2024 16:49:27 +0100 Subject: [PATCH 62/62] Changes to incorporate prognostic B22 --- src/2d/parcels/parcel_bc.f90 | 2 +- src/2d/parcels/parcel_container.f90 | 4 ++-- src/2d/parcels/parcel_diagnostics.f90 | 7 +++---- src/2d/parcels/parcel_ellipse.f90 | 21 +++++++-------------- src/2d/parcels/parcel_init.f90 | 10 ++++++---- src/2d/parcels/parcel_interpl.f90 | 9 ++++----- src/2d/parcels/parcel_merge.f90 | 13 +++++-------- src/2d/parcels/parcel_netcdf.f90 | 19 ++++++++++++++++++- src/2d/parcels/parcel_split.f90 | 3 ++- src/2d/stepper/ls_rk4.f90 | 16 +++++++++++++--- src/2d/stepper/rk4_utils.f90 | 11 +++++------ unit-tests/2d/test_ellipse_orientation.f90 | 7 +++++-- unit-tests/2d/test_ellipse_reflection.f90 | 7 +++++-- 13 files changed, 76 insertions(+), 53 deletions(-) diff --git a/src/2d/parcels/parcel_bc.f90 b/src/2d/parcels/parcel_bc.f90 index a305c1aa9..40b737f62 100644 --- a/src/2d/parcels/parcel_bc.f90 +++ b/src/2d/parcels/parcel_bc.f90 @@ -23,7 +23,7 @@ end subroutine apply_periodic_bc ! @param[inout] position vector of parcel ! @param[inout] B matrix of parcel subroutine apply_reflective_bc(position, B) - double precision, intent(inout) :: position(2), B(2) + double precision, intent(inout) :: position(2), B(3) if (position(2) > upper(2)) then position(2) = two * upper(2) - position(2) diff --git a/src/2d/parcels/parcel_container.f90 b/src/2d/parcels/parcel_container.f90 index d2a17a3d5..0c698dd67 100644 --- a/src/2d/parcels/parcel_container.f90 +++ b/src/2d/parcels/parcel_container.f90 @@ -12,7 +12,7 @@ module parcel_container type parcel_container_type double precision, allocatable, dimension(:, :) :: & position, & - B ! B matrix entries; ordering B(:, 1) = B11, B(:, 2) = B12 + B ! B matrix entries; ordering B(1, :) = B11, B(2, :) = B12, B(3, :) = B33 double precision, allocatable, dimension(:) :: & volume, & @@ -88,7 +88,7 @@ subroutine parcel_alloc(num) allocate(parcels%position(2, num)) allocate(parcels%vorticity(num)) - allocate(parcels%B(2, num)) + allocate(parcels%B(3, num)) allocate(parcels%volume(num)) allocate(parcels%buoyancy(num)) #ifndef ENABLE_DRY_MODE diff --git a/src/2d/parcels/parcel_diagnostics.f90 b/src/2d/parcels/parcel_diagnostics.f90 index 31c9fcf98..460d74325 100644 --- a/src/2d/parcels/parcel_diagnostics.f90 +++ b/src/2d/parcels/parcel_diagnostics.f90 @@ -88,7 +88,7 @@ subroutine calculate_parcel_diagnostics(velocity) double precision :: velocity(:, :) integer :: n double precision :: b, z, vel(2), vol, zmin - double precision :: eval, lam, B22, lsum, l2sum, vsum, v2sum + double precision :: eval, lam, lsum, l2sum, vsum, v2sum call start_timer(parcel_stats_timer) @@ -121,7 +121,7 @@ subroutine calculate_parcel_diagnostics(velocity) std_vol = zero !$omp parallel default(shared) - !$omp do private(n, vel, vol, b, z, eval, lam, B22) & + !$omp do private(n, vel, vol, b, z, eval, lam) & !$omp& reduction(+: ke, ape, lsum, l2sum, vsum, v2sum, n_small, rms_zeta) & !$omp& reduction(-: pe) do n = 1, n_parcels @@ -141,8 +141,7 @@ subroutine calculate_parcel_diagnostics(velocity) ape = ape + ape_den(b, z) * vol endif - B22 = get_B22(parcels%B(1, n), parcels%B(2, n), vol) - eval = get_eigenvalue(parcels%B(1, n), parcels%B(2, n), B22) + eval = get_eigenvalue(parcels%B(1, n), parcels%B(2, n), parcels%B(3, n)) lam = get_aspect_ratio(eval, vol) lsum = lsum + lam diff --git a/src/2d/parcels/parcel_ellipse.f90 b/src/2d/parcels/parcel_ellipse.f90 index 015de2d7f..df9b85b88 100644 --- a/src/2d/parcels/parcel_ellipse.f90 +++ b/src/2d/parcels/parcel_ellipse.f90 @@ -55,17 +55,14 @@ function get_eigenvector(a2, B11, B12, B22) result(evec) end function get_eigenvector ! used in unit tests only - function get_angle(B11, B12, volume) result(angle) + function get_angle(B11, B12, B22) result(angle) double precision, intent(in) :: B11 double precision, intent(in) :: B12 - double precision, intent(in) :: volume - double precision :: B22 + double precision, intent(in) :: B22 double precision :: a2 double precision :: evec(2) double precision :: angle - B22 = get_B22(B11, B12, volume) - a2 = get_eigenvalue(B11, B12, B22) evec = get_eigenvector(a2, B11, B12, B22) @@ -119,24 +116,20 @@ elemental function get_aspect_ratio(a2, volume) result(lam) ! @param[in] volume of the parcel ! @param[in] B matrix elements of the parcel ! @returns the parcel support points - function get_ellipse_points(position, volume, B) result(points) + function get_ellipse_points(position, B) result(points) double precision, intent(in) :: position(2) - double precision, intent(in) :: volume - double precision, intent(in) :: B(2) ! B11, B12 - double precision :: B22 + double precision, intent(in) :: B(3) ! B11, B12, B22 double precision :: c double precision :: a2 double precision :: evec(2) double precision :: h(2) double precision :: points(2, 2) - B22 = get_B22(B(1), B(2), volume) - - a2 = get_eigenvalue(B(1), B(2), B22) + a2 = get_eigenvalue(B(1), B(2), B(3)) - c = dsqrt(dabs(two * a2 - B(1) - B22)) + c = dsqrt(dabs(two * a2 - B(1) - B(3))) - evec = get_eigenvector(a2, B(1), B(2), B22) + evec = get_eigenvector(a2, B(1), B(2), B(3)) h = f12 * c * evec diff --git a/src/2d/parcels/parcel_init.f90 b/src/2d/parcels/parcel_init.f90 index 2e1c04ca4..72fd49af0 100644 --- a/src/2d/parcels/parcel_init.f90 +++ b/src/2d/parcels/parcel_init.f90 @@ -5,7 +5,7 @@ module parcel_init use options, only : parcel, output, verbose, field_tol use constants, only : zero, two, one, f12 use parcel_container, only : parcels, n_parcels - use parcel_ellipse, only : get_ab, get_B22, get_eigenvalue + use parcel_ellipse, only : get_ab, get_eigenvalue, get_B22 use parcel_split, only : split_ellipses use parcel_interpl, only : bilinear, ngp use parameters, only : dx, vcell, ncell, & @@ -83,6 +83,9 @@ subroutine init_parcels(fname, tol) ! B12 parcels%B(2, n) = zero + + ! B22 + parcels%B(3, n) = get_B22(parcels%B(1, n), zero, parcels%volume(1)) enddo !$omp end do !$omp end parallel @@ -145,13 +148,12 @@ end subroutine init_regular_positions subroutine init_refine(lam) double precision, intent(inout) :: lam - double precision :: B22, a2 + double precision :: a2 ! do refining by splitting do while (lam >= parcel%lambda_max) call split_ellipses(parcels, parcel%lambda_max) - B22 = get_B22(parcels%B(1, 1), zero, parcels%volume(1)) - a2 = get_eigenvalue(parcels%B(1, 1), zero, B22) + a2 = get_eigenvalue(parcels%B(1, 1), parcels%B(2, 1), parcels%B(3, 1)) lam = a2 / get_ab(parcels%volume(1)) end do end subroutine init_refine diff --git a/src/2d/parcels/parcel_interpl.f90 b/src/2d/parcels/parcel_interpl.f90 index 1f5d7783a..6abb54ab1 100644 --- a/src/2d/parcels/parcel_interpl.f90 +++ b/src/2d/parcels/parcel_interpl.f90 @@ -50,7 +50,7 @@ subroutine vol2grid pvol = parcels%volume(n) points = get_ellipse_points(parcels%position(:, n), & - pvol, parcels%B(:, n)) + parcels%B(:, n)) ! we have 2 points per ellipse @@ -85,7 +85,7 @@ end subroutine vol2grid #ifndef NDEBUG ! Interpolate the parcel volume to the grid to check symmetry subroutine vol2grid_symmetry_error - double precision :: points(2, 2), V, B(2), pos(2) + double precision :: points(2, 2), V, B(3), pos(2) integer :: n, p, l, m double precision :: pvol @@ -107,7 +107,7 @@ subroutine vol2grid_symmetry_error B(2) = dble(m) * B(2) - points = get_ellipse_points(pos, V, B) + points = get_ellipse_points(pos, B) ! we have 2 points per ellipse do p = 1, 2 @@ -180,7 +180,7 @@ subroutine par2grid btot = parcels%buoyancy(n) #endif points = get_ellipse_points(parcels%position(:, n), & - pvol, parcels%B(:, n)) + parcels%B(:, n)) call get_index(parcels%position(:, n), i, j) i = mod(i + nx, nx) @@ -334,7 +334,6 @@ subroutine grid2par(vel, vor, vgrad, add) do n = 1, n_parcels points = get_ellipse_points(parcels%position(:, n), & - parcels%volume(n), & parcels%B(:, n)) ! we have 2 points per ellipse diff --git a/src/2d/parcels/parcel_merge.f90 b/src/2d/parcels/parcel_merge.f90 index 8459aafb7..9f8a5f4ad 100644 --- a/src/2d/parcels/parcel_merge.f90 +++ b/src/2d/parcels/parcel_merge.f90 @@ -9,7 +9,7 @@ module parcel_merge , n_parcels & , parcel_replace & , get_delx - use parcel_ellipse, only : get_B22, get_ab + use parcel_ellipse, only : get_ab use options, only : parcel, verbose use parcel_bc use timer, only : start_timer, stop_timer @@ -80,7 +80,7 @@ subroutine do_group_merge(parcels, isma, iclo, n_merge, B11m, B12m, B22m, vm) integer :: m, ic, is, l, n integer :: loca(n_parcels) double precision :: x0(n_merge) - double precision :: posm(2, n_merge), delx, vmerge, dely, B22, mu + double precision :: posm(2, n_merge), delx, vmerge, dely, mu double precision :: buoym(n_merge), vortm(n_merge) #ifndef ENABLE_DRY_MODE double precision :: hum(n_merge) @@ -186,15 +186,13 @@ subroutine do_group_merge(parcels, isma, iclo, n_merge, B11m, B12m, B22m, vm) vmerge = one / vm(l) - B22 = get_B22(parcels%B(1, ic), parcels%B(2, ic), parcels%volume(ic)) - delx = get_delx(parcels%position(1, ic), posm(1, l)) dely = parcels%position(2, ic) - posm(2, l) mu = parcels%volume(ic) * vmerge B11m(l) = mu * (four * delx ** 2 + parcels%B(1, ic)) B12m(l) = mu * (four * delx * dely + parcels%B(2, ic)) - B22m(l) = mu * (four * dely ** 2 + B22) + B22m(l) = mu * (four * dely ** 2 + parcels%B(3, ic)) parcels%volume(ic) = vm(l) parcels%position(1, ic) = posm(1, l) @@ -216,14 +214,12 @@ subroutine do_group_merge(parcels, isma, iclo, n_merge, B11m, B12m, B22m, vm) delx = get_delx(parcels%position(1, is), posm(1, n)) dely = parcels%position(2, is) - posm(2, n) - B22 = get_B22(parcels%B(1, is), parcels%B(2, is), parcels%volume(is)) - ! volume fraction A_{is} / A mu = vmerge * parcels%volume(is) B11m(n) = B11m(n) + mu * (four * delx ** 2 + parcels%B(1, is)) B12m(n) = B12m(n) + mu * (four * delx * dely + parcels%B(2, is)) - B22m(n) = B22m(n) + mu * (four * dely ** 2 + B22) + B22m(n) = B22m(n) + mu * (four * dely ** 2 + parcels%B(3, is)) enddo end subroutine do_group_merge @@ -266,6 +262,7 @@ subroutine geometric_merge(parcels, isma, iclo, n_merge) parcels%B(1, ic) = B11(l) * factor parcels%B(2, ic) = B12(l) * factor + parcels%B(3, ic) = B22(l) * factor call apply_periodic_bc(parcels%position(:, ic)) endif diff --git a/src/2d/parcels/parcel_netcdf.f90 b/src/2d/parcels/parcel_netcdf.f90 index 9b86c1df9..cae7ca7df 100644 --- a/src/2d/parcels/parcel_netcdf.f90 +++ b/src/2d/parcels/parcel_netcdf.f90 @@ -21,7 +21,7 @@ module parcel_netcdf integer :: ncid integer :: npar_dim_id, vol_id, buo_id, & x_pos_id, z_pos_id, vor_id, & - b11_id, b12_id, & + b11_id, b12_id, b22_id, & t_axis_id, t_dim_id double precision :: restart_time @@ -133,6 +133,15 @@ subroutine create_netcdf_parcel_file(basename, overwrite, l_restart) dimids=dimids, & varid=b12_id) + call define_netcdf_dataset(ncid=ncid, & + name='B22', & + long_name='B22 element of shape matrix', & + std_name='', & + unit='m^2', & + dtype=NF90_DOUBLE, & + dimids=dimids, & + varid=b22_id) + call define_netcdf_dataset(ncid=ncid, & name='volume', & long_name='parcel volume', & @@ -205,6 +214,7 @@ subroutine write_netcdf_parcels(t) call write_netcdf_dataset(ncid, b11_id, parcels%B(1, 1:n_parcels), start, cnt) call write_netcdf_dataset(ncid, b12_id, parcels%B(2, 1:n_parcels), start, cnt) + call write_netcdf_dataset(ncid, b22_id, parcels%B(3, 1:n_parcels), start, cnt) call write_netcdf_dataset(ncid, vol_id, parcels%volume(1:n_parcels), start, cnt) @@ -263,6 +273,13 @@ subroutine read_netcdf_parcels(fname) stop endif + if (has_dataset(ncid, 'B22')) then + call read_netcdf_dataset(ncid, 'B22', parcels%B(3, 1:n_parcels), start, cnt) + else + print *, "The parcel shape component B22 must be present! Exiting." + stop + endif + if (has_dataset(ncid, 'x_position')) then call read_netcdf_dataset(ncid, 'x_position', & parcels%position(1, 1:n_parcels), start, cnt) diff --git a/src/2d/parcels/parcel_split.f90 b/src/2d/parcels/parcel_split.f90 index 59e7d2967..f8a738efa 100644 --- a/src/2d/parcels/parcel_split.f90 +++ b/src/2d/parcels/parcel_split.f90 @@ -45,7 +45,7 @@ subroutine split_ellipses(parcels, threshold) B11 = parcels%B(1, n) B12 = parcels%B(2, n) V = parcels%volume(n) - B22 = get_B22(B11, B12, V) + B22 = parcels%B(3, n) a2 = get_eigenvalue(B11, B12, B22) @@ -64,6 +64,7 @@ subroutine split_ellipses(parcels, threshold) parcels%B(1, n) = B11 - f34 * a2 * evec(1) ** 2 parcels%B(2, n) = B12 - f34 * a2 * (evec(1) * evec(2)) + parcels%B(3, n) = B22 - f34 * a2 * evec(2) ** 2 h = f14 * dsqrt(three * a2) parcels%volume(n) = f12 * V diff --git a/src/2d/stepper/ls_rk4.f90 b/src/2d/stepper/ls_rk4.f90 index bd072f1e9..4468a90b9 100644 --- a/src/2d/stepper/ls_rk4.f90 +++ b/src/2d/stepper/ls_rk4.f90 @@ -8,6 +8,7 @@ module ls_rk4 use parcel_bc use rk4_utils, only: evolve_ellipsoid, get_time_step use utils, only : write_step + use parcel_ellipse, only : get_ab use parcel_interpl, only : par2grid, grid2par, grid2par_add use fields, only : velgradg, velog, vortg, vtend, tbuoyg use tri_inversion, only : vor2vel, vorticity_tendency @@ -122,7 +123,7 @@ end subroutine ls_rk4_step subroutine ls_rk4_substep(dt, step) double precision, intent(in) :: dt integer, intent(in) :: step - double precision :: ca, cb + double precision :: ca, cb, factor integer :: n ca = cas(step) @@ -152,14 +153,23 @@ subroutine ls_rk4_substep(dt, step) + cb * dt * delta_pos(:, n) parcels%vorticity(n) = parcels%vorticity(n) + cb * dt * delta_vor(n) - call evolve_ellipsoid(parcels%B(:, n), int_strain(:, n), parcels%volume(n), cb * dt) + call evolve_ellipsoid(parcels%B(:, n), int_strain(:, n), cb * dt) enddo !$omp end parallel do call stop_timer(rk4_timer) if (step == 5) then - return + ! Correct parcel volume + !$omp parallel do default(shared) private(n, factor) + do n = 1, n_parcels + factor = get_ab(parcels%volume(n)) / dsqrt(parcels%B(1, n) * parcels%B(3, n) - parcels%B(2, n) ** 2) + parcels%B(1, n) = parcels%B(1, n) * factor + parcels%B(2, n) = parcels%B(2, n) * factor + parcels%B(3, n) = parcels%B(3, n) * factor + end do + !$omp end parallel do + return endif call start_timer(rk4_timer) diff --git a/src/2d/stepper/rk4_utils.f90 b/src/2d/stepper/rk4_utils.f90 index 971d2f2e3..ee03712ff 100644 --- a/src/2d/stepper/rk4_utils.f90 +++ b/src/2d/stepper/rk4_utils.f90 @@ -1,5 +1,4 @@ module rk4_utils - use parcel_ellipse, only : get_B22 use fields, only : velgradg, tbuoyg, vtend use constants, only : zero, one, two, f12, f14, f23, c_matexp_x1, c_matexp_x2, c_matexp_x4, & c_matexp_x5, c_matexp_x6, c_matexp_x7, c_matexp_y2 @@ -13,10 +12,9 @@ module rk4_utils contains - subroutine evolve_ellipsoid(B, S, volume, dt_sub) - double precision, intent(inout) :: B(2) + subroutine evolve_ellipsoid(B, S, dt_sub) + double precision, intent(inout) :: B(3) double precision, intent(in) :: S(4) - double precision, intent(in) :: volume double precision, intent(in) :: dt_sub double precision :: Bmat(2,2) double precision :: Smat(2,2) @@ -28,8 +26,8 @@ subroutine evolve_ellipsoid(B, S, volume, dt_sub) Bmat(1, 1) = B(1) ! B11 Bmat(1, 2) = B(2) ! B12 - Bmat(2, 1) = B(1) ! B21 - Bmat(2, 2) = get_B22(B(1), B(2), volume) + Bmat(2, 1) = B(2) ! B21 + Bmat(2, 2) = B(3) ! B22 Smat(1, 1) = S(1) ! S11 Smat(1, 2) = S(2) ! S12 @@ -61,6 +59,7 @@ subroutine evolve_ellipsoid(B, S, volume, dt_sub) B(1) = Bmat(1, 1) B(2) = Bmat(1, 2) + B(3) = Bmat(2, 2) end subroutine evolve_ellipsoid diff --git a/unit-tests/2d/test_ellipse_orientation.f90 b/unit-tests/2d/test_ellipse_orientation.f90 index 4fd3f123d..5b8cf6cf7 100644 --- a/unit-tests/2d/test_ellipse_orientation.f90 +++ b/unit-tests/2d/test_ellipse_orientation.f90 @@ -15,7 +15,7 @@ program test_ellipse_orientation double precision :: extent(2) = (/0.2d0, 0.2d0/) integer :: iter integer :: grid(2) = (/2, 2/) - double precision :: angle, B11, B12, V, a2, b2 + double precision :: angle, B11, B12, B22, V, a2, b2 double precision, parameter :: lam = three logical :: passed = .true. @@ -41,6 +41,9 @@ program test_ellipse_orientation parcels%B(1, 1) = B11 parcels%B(2, 1) = B12 + + B22 = get_B22(B11, B12, V) + parcels%B(3, 1) = B22 ! get_angle computes the angle in the first and fourth quadrant, i.e., ! -pi/2 <= get_angle <= pi/2 @@ -50,7 +53,7 @@ program test_ellipse_orientation angle = angle - two * pi endif - passed = (passed .and. abs(angle - get_angle(B11, B12, V)) < dble(1.0e-13)) + passed = (passed .and. abs(angle - get_angle(B11, B12, B22)) < dble(1.0e-13)) enddo diff --git a/unit-tests/2d/test_ellipse_reflection.f90 b/unit-tests/2d/test_ellipse_reflection.f90 index d35d56a99..4af4d8e6e 100644 --- a/unit-tests/2d/test_ellipse_reflection.f90 +++ b/unit-tests/2d/test_ellipse_reflection.f90 @@ -10,12 +10,12 @@ program test_ellipse_reflection use constants, only : pi, zero, one, two, three, four, f12, f14 use parcel_container use parcel_bc, only : apply_reflective_bc - use parcel_ellipse, only : get_ab, get_angle + use parcel_ellipse, only : get_ab, get_angle, get_B22 use parameters, only : update_parameters, lower, dx, vcell, extent, nx, nz implicit none integer :: iter - double precision :: angle, ab, B11, B12, error, a2, b2 + double precision :: angle, ab, B11, B12, B22, error, a2, b2 nx = 10 nz = 10 @@ -48,6 +48,9 @@ program test_ellipse_reflection parcels%B(1, 1) = B11 parcels%B(2, 1) = B12 + B22 = get_B22(B11, B12, V) + parcels%B(3, 1) = B22 + call apply_reflective_bc(parcels%position(:, 1), parcels%B(:, 1)) call check_result('lower')