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 ac171e276..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) @@ -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,10 +333,7 @@ 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)) ! 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 66a50855d..4468a90b9 100644 --- a/src/2d/stepper/ls_rk4.f90 +++ b/src/2d/stepper/ls_rk4.f90 @@ -6,8 +6,9 @@ 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_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 @@ -23,7 +24,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 +54,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 +64,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 +91,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) @@ -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) @@ -131,29 +132,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,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) - parcels%B(:, n) = parcels%B(:, n) + cb * dt * delta_b(:, n) + 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) @@ -181,7 +178,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..ee03712ff 100644 --- a/src/2d/stepper/rk4_utils.f90 +++ b/src/2d/stepper/rk4_utils.f90 @@ -1,34 +1,67 @@ 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, 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 #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, dt_sub) + double precision, intent(inout) :: B(3) 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) + 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 + Bmat(2, 1) = B(2) ! B21 + Bmat(2, 2) = B(3) ! B22 + + Smat(1, 1) = S(1) ! S11 + Smat(1, 2) = S(2) ! S12 + Smat(2, 1) = S(3) ! S21 + Smat(2, 2) = S(4) ! S33 + + ! 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 + + ! 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, 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))) + + B(1) = Bmat(1, 1) + B(2) = Bmat(1, 2) + B(3) = Bmat(2, 2) + + end subroutine evolve_ellipsoid ! Estimate a suitable time step based on the velocity strain ! and buoyancy gradient. 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) diff --git a/src/3d/parcels/parcel_container.f90 b/src/3d/parcels/parcel_container.f90 index bf273fd3b..8e54ea83d 100644 --- a/src/3d/parcels/parcel_container.f90 +++ b/src/3d/parcels/parcel_container.f90 @@ -37,6 +37,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 @@ -57,6 +58,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_DUDZ, & ! RK4 variable du/dz @@ -75,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(1, :) = B11, B(2, :) = B12, B(3, :) = B13 + ! B(4, :) = B22, B(5, :) = B23, B(6, :) = B33 double precision, allocatable, dimension(:) :: & volume, & @@ -92,8 +94,7 @@ module parcel_container double precision, allocatable, dimension(:, :) :: & delta_pos, & ! velocity delta_vor, & ! vorticity tendency - strain, & - delta_b ! B-matrix tendency + int_strain #ifdef ENABLE_LABELS integer(kind=8), allocatable, dimension(:) :: & @@ -124,8 +125,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 @@ -134,7 +136,7 @@ subroutine set_buffer_indices #endif #ifdef ENABLE_LABELS IDX_LABEL = i - IDX_DILUTION = i +1 + IDX_DILUTION = i + 1 i = i + 2 #endif @@ -145,21 +147,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_DUDX = i + 11 - IDX_RK4_DUDY = i + 12 - IDX_RK4_DUDZ = i + 13 - IDX_RK4_DVDX = i + 14 - IDX_RK4_DVDY = i + 15 - IDX_RK4_DVDZ = i + 16 - IDX_RK4_DWDX = i + 17 - IDX_RK4_DWDY = i + 18 - - i = i + 19 + 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) @@ -262,12 +259,12 @@ 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 #ifdef ENABLE_LABELS parcels%label(n) = parcels%label(m) @@ -278,10 +275,9 @@ subroutine parcel_replace(n, 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%delta_b(:, n) = parcels%delta_b(:, 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 @@ -319,8 +315,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%delta_b, new_size, n_parcels) + call resize_array(parcels%int_strain, new_size, n_parcels) call stop_timer(resize_timer) @@ -337,7 +332,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 @@ -352,8 +347,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%delta_b(5, num)) + allocate(parcels%int_strain(8, num)) end subroutine parcel_alloc @@ -386,8 +380,7 @@ subroutine parcel_dealloc ! LS-RK4 variables deallocate(parcels%delta_pos) deallocate(parcels%delta_vor) - deallocate(parcels%strain) - deallocate(parcels%delta_b) + deallocate(parcels%int_strain) end subroutine parcel_dealloc @@ -400,7 +393,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 @@ -413,8 +406,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_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 @@ -428,7 +420,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 @@ -441,8 +433,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%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_damping.f90 b/src/3d/parcels/parcel_damping.f90 index 1e0b96779..064fa4d42 100644 --- a/src/3d/parcels/parcel_damping.f90 +++ b/src/3d/parcels/parcel_damping.f90 @@ -125,7 +125,7 @@ subroutine perturbation_damping(dt, l_reuse) pvol = parcels%volume(n) #ifndef ENABLE_P2G_1POINT points = get_ellipsoid_points(parcels%position(:, n), & - pvol, parcels%B(:, n), n, l_reuse) + parcels%B(:, n), n, l_reuse) #else points(:, 1) = parcels%position(:, n) #endif diff --git a/src/3d/parcels/parcel_diagnostics.f90 b/src/3d/parcels/parcel_diagnostics.f90 index 3879d4938..7aa8b884a 100644 --- a/src/3d/parcels/parcel_diagnostics.f90 +++ b/src/3d/parcels/parcel_diagnostics.f90 @@ -87,7 +87,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 @@ -106,7 +106,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 diff --git a/src/3d/parcels/parcel_ellipsoid.f90 b/src/3d/parcels/parcel_ellipsoid.f90 index 03ec2f208..afdd7c250 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 @@ -129,12 +130,13 @@ 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, 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,22 +147,22 @@ 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 !:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: ! 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, volume) result(D) - double precision, intent(in) :: B(5) - double precision, intent(in) :: volume + function get_eigenvalues(B) result(D) + double precision, intent(in) :: B(6) + ! 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) @@ -175,15 +177,13 @@ 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 - - B33 = get_B33(B, volume) + function get_determinant(B) result(detB) + double precision, intent(in) :: B(6) + ! double precision, intent(in) :: volume + double precision :: detB - 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 +195,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 +224,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) - double precision, intent(in) :: volume + 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) @@ -301,10 +301,10 @@ 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) :: B(5) ! B11, B12, B13, B22, B23 + ! double precision, intent(in) :: volume + 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) @@ -317,7 +317,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) @@ -327,7 +327,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) @@ -351,13 +351,13 @@ end function get_ellipsoid_points !:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: - function get_angles(B, volume) result(angles) - double precision, intent(in) :: B(5) - double precision, intent(in) :: volume + function get_angles(B) result(angles) + double precision, intent(in) :: B(6) + ! 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_init.f90 b/src/3d/parcels/parcel_init.f90 index 8da1bee24..5769cacdf 100644 --- a/src/3d/parcels/parcel_init.f90 +++ b/src/3d/parcels/parcel_init.f90 @@ -88,6 +88,9 @@ subroutine parcel_default ! B22 parcels%B(4, n) = l23 + + ! B33 + parcels%B(6, n) = l23 enddo !$omp end do !$omp end parallel @@ -168,7 +171,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 diff --git a/src/3d/parcels/parcel_interpl.f90 b/src/3d/parcels/parcel_interpl.f90 index 15a27c282..30f292e4a 100644 --- a/src/3d/parcels/parcel_interpl.f90 +++ b/src/3d/parcels/parcel_interpl.f90 @@ -107,7 +107,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 @@ -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,14 @@ subroutine par2grid(l_reuse) call start_timer(par2grid_timer) + ! This is only here to allow debug compilation + ! with a warning for unused variables +#if defined (ENABLE_P2G_1POINT) && !defined (NDEBUG) + if(present(l_reuse)) then + l_reuse=.false. + endif +#endif + vortg = zero volg = zero nparg = zero @@ -200,7 +209,7 @@ subroutine par2grid(l_reuse) #ifndef ENABLE_P2G_1POINT points = get_ellipsoid_points(parcels%position(:, n), & - pvol, parcels%B(:, n), n, l_reuse) + parcels%B(:, n), n, l_reuse) #else points(:, 1) = parcels%position(:, n) #endif @@ -423,6 +432,7 @@ subroutine grid2par(add) do n = 1, n_parcels parcels%delta_pos(:, n) = zero parcels%delta_vor(:, n) = zero + parcels%int_strain(:, n) = zero enddo !$omp end do !$omp end parallel @@ -433,6 +443,7 @@ subroutine grid2par(add) do n = 1, n_parcels parcels%delta_pos(:, n) = zero parcels%delta_vor(:, n) = zero + parcels%int_strain(:, n) = zero enddo !$omp end do !$omp end parallel @@ -442,11 +453,9 @@ subroutine grid2par(add) !$omp do private(n, l, p, points, is, js, ks, weights) do n = 1, n_parcels - parcels%strain(:, n) = zero - #ifndef ENABLE_G2P_1POINT points = get_ellipsoid_points(parcels%position(:, n), & - parcels%volume(n), parcels%B(:, n), n) + parcels%B(:, n), n) #else points(:, 1) = parcels%position(:, n) #endif @@ -460,8 +469,8 @@ subroutine grid2par(add) + point_weight_g2p * sum(weights * velog(ks:ks+1, js:js+1, is:is+1, l)) enddo do l = 1,8 - parcels%strain(l, n) = parcels%strain(l, n) & - + point_weight_g2p * sum(weights * velgradg(ks:ks+1, js:js+1, is:is+1, l)) + 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)) enddo do l = 1,3 parcels%delta_vor(l, n) = parcels%delta_vor(l, n) & diff --git a/src/3d/parcels/parcel_merge.f90 b/src/3d/parcels/parcel_merge.f90 index 08cd23c9f..49adc2609 100644 --- a/src/3d/parcels/parcel_merge.f90 +++ b/src/3d/parcels/parcel_merge.f90 @@ -13,7 +13,7 @@ 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 options, only : parcel use datatypes, only : int64 #if defined (ENABLE_VERBOSE) && !defined (NDEBUG) @@ -127,7 +127,7 @@ subroutine do_group_merge(isma, iclo, n_merge, Bm, vm) integer :: m, ic, is, l, n 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) @@ -268,8 +268,6 @@ subroutine do_group_merge(isma, iclo, n_merge, Bm, vm) vmerge = one / vm(l) - 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)) delz = parcels%position(3, ic) - posm(3, l) @@ -282,7 +280,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) @@ -310,8 +308,6 @@ 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)) - ! volume fraction V_{is} / V mu = vmerge * parcels%volume(is) @@ -320,7 +316,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(6, is)) enddo end subroutine do_group_merge @@ -362,7 +358,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_netcdf.f90 b/src/3d/parcels/parcel_netcdf.f90 index 63526bc4f..e8c025b5e 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) @@ -606,6 +608,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) @@ -722,6 +731,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) @@ -768,6 +778,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 @@ -790,7 +801,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) @@ -896,6 +907,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='', & diff --git a/src/3d/parcels/parcel_split.f90 b/src/3d/parcels/parcel_split.f90 index af2a6186b..fe481412c 100644 --- a/src/3d/parcels/parcel_split.f90 +++ b/src/3d/parcels/parcel_split.f90 @@ -39,7 +39,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 @@ -57,12 +57,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) @@ -118,7 +118,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 @@ -130,6 +130,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 diff --git a/src/3d/stepper/ls_rk.f90 b/src/3d/stepper/ls_rk.f90 index 0bb5529df..fc8dd7df6 100644 --- a/src/3d/stepper/ls_rk.f90 +++ b/src/3d/stepper/ls_rk.f90 @@ -3,12 +3,14 @@ ! (see https://doi.org/10.5194/gmd-10-3145-2017) ! ============================================================================= module ls_rk + use constants, only : f13, one use options, only : time use dimensions, only : I_Z use parcel_container use parcel_bc use parcel_mpi, only : parcel_communicate - use rk_utils, only: get_dBdt, get_time_step + use parcel_ellipsoid, only : get_abc + use rk_utils, only: get_time_step, evolve_ellipsoid use utils, only : write_step use parcel_interpl, only : par2grid, grid2par use parcel_damping, only : parcel_damp @@ -24,7 +26,6 @@ module ls_rk private integer, parameter :: dp=kind(zero) ! double precision - integer :: rk_timer ! fourth order RK coefficients: @@ -122,15 +123,16 @@ 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) 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 @@ -143,42 +145,18 @@ 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) cb = cbptr(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), & - parcels%volume(n)) - enddo - !$omp end parallel do - - call stop_timer(rk_timer) - else + if (step > 1) then call vor2vel call vorticity_tendency call grid2par(add=.true.) - - 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), & - parcels%volume(n)) - enddo - !$omp end parallel do - - call stop_timer(rk_timer) endif call start_timer(rk_timer) @@ -190,25 +168,38 @@ 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) + + call evolve_ellipsoid(parcels%B(:, n), parcels%int_strain(:, n), cb * dt) + 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 - return - endif - call start_timer(rk_timer) + call stop_timer(rk_timer) + + return + endif !$omp parallel do default(shared) private(n) 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%int_strain(:, n) = ca * parcels%int_strain(:, n) enddo !$omp end parallel do diff --git a/src/3d/stepper/rk_utils.f90 b/src/3d/stepper/rk_utils.f90 index e3901a547..efed42c3c 100644 --- a/src/3d/stepper/rk_utils.f90 +++ b/src/3d/stepper/rk_utils.f90 @@ -1,9 +1,10 @@ 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_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, 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 @@ -17,54 +18,67 @@ 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 - ! 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, volume) result(Bout) - double precision, intent(in) :: Bin(I_B23) + subroutine evolve_ellipsoid(B, S, dt_sub) + double precision, intent(inout) :: B(I_B33) double precision, intent(in) :: S(8) - double precision, intent(in) :: volume - double precision :: Bout(5), 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) * 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) * B33 ! + dv/dz * B33 - end function get_dBdt + 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 + Bmat(1, 2) = B(I_B12) ! B12 + Bmat(1, 3) = B(I_B13) ! B13 + Bmat(2, 1) = Bmat(1, 2) ! B21 + Bmat(2, 2) = B(I_B22) ! B22 + Bmat(2, 3) = B(I_B23) ! B23 + 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 + 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) = -(Smat(1, 1) + Smat(2, 2)) ! S33 + + ! 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, 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))) + + 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 + ! Calculate velocity strain ! @param[in] velocity gradient tensor at grid point diff --git a/src/utils/constants.f90 b/src/utils/constants.f90 index 5cfd74699..e92256da2 100644 --- a/src/utils/constants.f90 +++ b/src/utils/constants.f90 @@ -45,4 +45,19 @@ module constants double precision, parameter :: rad2deg = 180.0d0 * fpi double precision, parameter :: deg2rad = one / rad2deg + ! 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 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') diff --git a/unit-tests/3d/test_mpi_gradient_correction_3d.f90 b/unit-tests/3d/test_mpi_gradient_correction_3d.f90 index b8286e128..7be156473 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,11 @@ program test_mpi_gradient_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_grid2par.f90 b/unit-tests/3d/test_mpi_grid2par.f90 index 6f72c1b6f..2b1572a20 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. @@ -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 @@ -98,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) 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) 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_diagnostics.f90 b/unit-tests/3d/test_mpi_parcel_diagnostics.f90 index 3496b9a05..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 @@ -66,6 +67,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 f1f41f63d..436eb0a74 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 @@ -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..bcfedd62e 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(6, 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..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) @@ -141,7 +142,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 +154,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..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 @@ -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. @@ -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)