m_perturbation.fpp.f90

Go to the documentation of this file.

# 1 "/home/runner/work/MFC/MFC/src/pre_process/m_perturbation.fpp"
!>
!! @file
!! @brief Contains module m_perturbation
 
!> @brief Perturbs initial mean flow fields with random noise, mixing-layer instabilities, or simplex noise
module m_perturbation
 
    use m_derived_types
    use m_global_parameters
    use m_mpi_proxy
    use m_boundary_common
    use m_helper
    use m_simplex_noise
    use ieee_arithmetic
 
    implicit none
 
    real(wp), allocatable, dimension(:,:,:,:) :: q_prim_temp
 
contains
 
    !> Allocate the temporary primitive variable array used by elliptic smoothing.
    impure subroutine s_initialize_perturbation_module()
 
        if (elliptic_smoothing) then
            allocate (q_prim_temp(0:m,0:n,0:p,1:sys_size))
        end if
 
    end subroutine s_initialize_perturbation_module
 
    !> Randomly perturb partial density fields at the interface of a spherical volume fraction region.
    impure subroutine s_perturb_sphere(q_prim_vf)
 
        type(scalar_field), dimension(sys_size), intent(inout) :: q_prim_vf
        integer                                                :: i, j, k, l
        real(wp)                                               :: perturb_alpha
        real(wp)                                               :: rand_real
 
        call random_seed()
 
        do k = 0, p
            do j = 0, n
                do i = 0, m
                    call random_number(rand_real)
 
                    perturb_alpha = q_prim_vf(eqn_idx%E + perturb_sph_fluid)%sf(i, j, k)
 
                    ! Perturb partial density fields to match perturbed volume fraction fields when the volume fraction is not near
                    ! 0 or 1
                    if ((.not. f_approx_equal(perturb_alpha, 0._wp)) .and. (.not. f_approx_equal(perturb_alpha, 1._wp))) then
                        do l = 1, num_fluids
                            q_prim_vf(l)%sf(i, j, k) = q_prim_vf(eqn_idx%E + l)%sf(i, j, k)*fluid_rho(l)
                        end do
                    end if
                end do
            end do
        end do
 
    end subroutine s_perturb_sphere
 
    !> Add random noise to the velocity and void fraction of the surrounding flow field.
    impure subroutine s_perturb_surrounding_flow(q_prim_vf)
 
        type(scalar_field), dimension(sys_size), intent(inout) :: q_prim_vf
        integer                                                :: i, j, k
        real(wp)                                               :: rand_real
 
        call random_seed()
 
        do k = 0, p
            do j = 0, n
                do i = 0, m
                    call random_number(rand_real)
                    rand_real = rand_real*perturb_flow_mag
                    q_prim_vf(eqn_idx%mom%end)%sf(i, j, k) = rand_real*q_prim_vf(eqn_idx%mom%beg)%sf(i, j, k)
                    q_prim_vf(eqn_idx%mom%beg)%sf(i, j, k) = (1._wp + rand_real)*q_prim_vf(eqn_idx%mom%beg)%sf(i, j, k)
                    if (bubbles_euler) then
                        q_prim_vf(eqn_idx%alf)%sf(i, j, k) = (1._wp + rand_real)*q_prim_vf(eqn_idx%alf)%sf(i, j, k)
                    end if
                end do
            end do
        end do
 
    end subroutine s_perturb_surrounding_flow
 
    !> Iteratively smooth all primitive variable fields using a discrete elliptic (Laplacian) filter.
    impure subroutine s_elliptic_smoothing(q_prim_vf, bc_type, q_T_sf)
 
        type(scalar_field), dimension(sys_size), intent(inout)     :: q_prim_vf
        type(integer_field), dimension(1:num_dims,1:2), intent(in) :: bc_type
        type(scalar_field), optional, intent(inout)                :: q_t_sf
        integer                                                    :: i, j, k, l, q
 
        do q = 1, elliptic_smoothing_iters
            ! Communication of buffer regions and apply boundary conditions
            call s_populate_variables_buffers(bc_type, q_prim_vf, pb%sf, mv%sf, q_t_sf)
 
            ! Perform smoothing and store in temp array
            if (n == 0) then
                do j = 0, m
                    do i = 1, sys_size
                        q_prim_temp(j, 0, 0, i) = (1._wp/4._wp)*(q_prim_vf(i)%sf(j + 1, 0, 0) + q_prim_vf(i)%sf(j - 1, 0, &
                                    & 0) + 2._wp*q_prim_vf(i)%sf(j, 0, 0))
                    end do
                end do
            else if (p == 0) then
                do k = 0, n
                    do j = 0, m
                        do i = 1, sys_size
                            q_prim_temp(j, k, 0, i) = (1._wp/8._wp)*(q_prim_vf(i)%sf(j + 1, k, 0) + q_prim_vf(i)%sf(j - 1, k, &
                                        & 0) + q_prim_vf(i)%sf(j, k + 1, 0) + q_prim_vf(i)%sf(j, k - 1, &
                                        & 0) + 4._wp*q_prim_vf(i)%sf(j, k, 0))
                        end do
                    end do
                end do
            else
                do l = 0, p
                    do k = 0, n
                        do j = 0, m
                            do i = 1, sys_size
                                q_prim_temp(j, k, l, i) = (1._wp/12._wp)*(q_prim_vf(i)%sf(j + 1, k, l) + q_prim_vf(i)%sf(j - 1, &
                                            & k, l) + q_prim_vf(i)%sf(j, k + 1, l) + q_prim_vf(i)%sf(j, k - 1, &
                                            & l) + q_prim_vf(i)%sf(j, k, l + 1) + q_prim_vf(i)%sf(j, k, &
                                            & l - 1) + 6._wp*q_prim_vf(i)%sf(j, k, l))
                            end do
                        end do
                    end do
                end do
            end if
 
            ! Copy smoothed data back to array of scalar fields
            do l = 0, p
                do k = 0, n
                    do j = 0, m
                        do i = 1, sys_size
                            q_prim_vf(i)%sf(j, k, l) = q_prim_temp(j, k, l, i)
                        end do
                    end do
                end do
            end do
        end do
 
    end subroutine s_elliptic_smoothing
 
    !> Perturb velocity and volume fraction fields using multi-octave simplex noise.
    subroutine s_perturb_simplex(q_prim_vf)
 
        type(scalar_field), dimension(sys_size), intent(inout) :: q_prim_vf
        real(wp)                                               :: mag, freq, scale, vel_rsm
        real(wp), dimension(:,:), allocatable                  :: ofs
        integer                                                :: nOffsets
        real(wp)                                               :: xl, yl, zl
        integer                                                :: i, j, k, l, q
 
        noffsets = max(num_dims, num_fluids)
 
        allocate (ofs(noffsets, num_dims))
 
        ! Store offsets
        do i = 1, num_dims
            do j = 1, num_dims
                ofs(j, i) = simplex_params%perturb_vel_offset(j, i)
            end do
        end do
 
        ! Perturb velocities
        do i = 1, num_dims
            if (simplex_params%perturb_vel(i)) then
                freq = simplex_params%perturb_vel_freq(i)
                scale = simplex_params%perturb_vel_scale(i)
                do l = 0, p
                    do k = 0, n
                        do j = 0, m
                            xl = freq*(x_cc(j) + ofs(i, 1))
                            yl = freq*(y_cc(k) + ofs(i, 2))
                            if (num_dims == 2) then
                                mag = f_simplex2d(xl, yl)
                            else if (num_dims == 3) then
                                zl = freq*(z_cc(l) + ofs(i, 3))
                                mag = f_simplex3d(xl, yl, zl)
                            end if
 
                            vel_rsm = 0._wp
                            do q = 1, num_dims
                                vel_rsm = vel_rsm + q_prim_vf(eqn_idx%mom%beg + q - 1)%sf(j, k, l)**2._wp
                            end do
                            vel_rsm = sqrt(vel_rsm)
 
                            q_prim_vf(eqn_idx%mom%beg + i - 1)%sf(j, k, l) = q_prim_vf(eqn_idx%mom%beg + i - 1)%sf(j, k, &
                                      & l) + vel_rsm*scale*mag
                        end do
                    end do
                end do
            end if
        end do
 
        ! Store offsets
        do i = 1, num_dims
            do j = 1, num_fluids
                ofs(j, i) = simplex_params%perturb_dens_offset(j, i)
            end do
        end do
 
        ! Perturb densities
        do i = 1, num_fluids
            if (simplex_params%perturb_dens(i)) then
                freq = simplex_params%perturb_dens_freq(i)
                scale = simplex_params%perturb_dens_scale(i)
                do l = 0, p
                    do k = 0, n
                        do j = 0, m
                            xl = freq*(x_cc(j) + ofs(i, 1))
                            yl = freq*(y_cc(k) + ofs(i, 2))
                            if (num_dims == 2) then
                                mag = f_simplex2d(xl, yl)
                            else if (num_dims == 3) then
                                zl = freq*(z_cc(l) + ofs(i, 3))
                                mag = f_simplex3d(xl, yl, zl)
                            end if
                            q_prim_vf(eqn_idx%cont%beg + i - 1)%sf(j, k, l) = q_prim_vf(eqn_idx%cont%beg + i - 1)%sf(j, k, &
                                      & l) + q_prim_vf(eqn_idx%cont%beg + i - 1)%sf(j, k, l)*scale*mag
                        end do
                    end do
                end do
            end if
        end do
 
        deallocate (ofs)
 
    end subroutine s_perturb_simplex
 
    !> Compute velocity perturbations for a temporal mixing layer with a hyperbolic tangent mean streamwise velocity profile, using
    !! an inverted version of the spectrum-based synthetic turbulence generation method proposed by Guo et al. (2023, JFM).
    subroutine s_perturb_mixlayer(q_prim_vf)
 
        type(scalar_field), dimension(sys_size), intent(inout) :: q_prim_vf
        real(wp), dimension(mixlayer_perturb_nk)               :: k, Ek
        real(wp), dimension(3, 3)                              :: Rij, Lmat
        real(wp), dimension(3)                                 :: velfluc, sig_tmp, sig, khat, xi
        real(wp)                                               :: dk, alpha, Eksum, q, uu0, phi
        integer                                                :: i, j, l, r, ierr
 
        dk = 1._wp/mixlayer_perturb_nk
 
        ! Compute prescribed energy spectra
        eksum = 0._wp
        do i = 1, mixlayer_perturb_nk
            k(i) = dk*i
            ek(i) = (k(i)/mixlayer_perturb_k0)**4._wp*exp(-2._wp*(k(i)/mixlayer_perturb_k0)**2._wp)
            eksum = eksum + ek(i)
        end do
 
        ! Main loop
        do r = 0, n
            ! Compute prescribed Reynolds stress tensor with about half magnitude of its self-similar value
            rij(:,:) = 0._wp
            uu0 = patch_icpp(1)%vel(1)**2._wp*(1._wp - tanh(y_cc(r)*mixlayer_vel_coef)**2._wp)
            rij(1, 1) = 0.05_wp*uu0
            rij(2, 2) = 0.03_wp*uu0
            rij(3, 3) = 0.03_wp*uu0
            rij(1, 2) = -0.02_wp*uu0
            rij(2, 1) = rij(1, 2)
 
            ! Cholesky decomposition for inhomogeneity and anisotropy
            lmat = 0._wp
            lmat(1, 1) = sqrt(rij(1, 1))
            if (abs(lmat(1, 1)) < sgm_eps) lmat(1, 1) = sgm_eps
            lmat(2, 1) = rij(2, 1)/lmat(1, 1)
            lmat(2, 2) = sqrt(rij(2, 2) - lmat(2, 1)**2._wp)
            if (abs(lmat(2, 2)) < sgm_eps) lmat(2, 2) = sgm_eps
            lmat(3, 1) = rij(3, 1)/lmat(1, 1)
            lmat(3, 2) = (rij(3, 2) - lmat(3, 1)*lmat(2, 1))/lmat(2, 2)
            lmat(3, 3) = sqrt(rij(3, 3) - lmat(3, 1)**2._wp - lmat(3, 2)**2._wp)
 
            ! Compute perturbation for each Fourier component
            do i = 1, mixlayer_perturb_nk
                ! Generate random numbers for unit wavevector khat, random unit vector xi, and random mode phase phi
                if (proc_rank == 0) then
                    call s_generate_random_perturbation(khat, xi, phi, i, y_cc(r))
                end if
 
#ifdef MFC_MPI
                call mpi_bcast(khat, 3, mpi_p, 0, mpi_comm_world, ierr)
                call mpi_bcast(xi, 3, mpi_p, 0, mpi_comm_world, ierr)
                call mpi_bcast(phi, 1, mpi_p, 0, mpi_comm_world, ierr)
#endif
 
                ! Compute mode direction by two-time cross product
                sig_tmp = f_cross(xi, khat)
                sig_tmp = sig_tmp/sqrt(sum(sig_tmp**2._wp))
                sig = f_cross(khat, sig_tmp)
 
                ! Compute perturbation for each grid
                do l = 0, p
                    do j = 0, m
                        q = sqrt(ek(i)/eksum)
                        alpha = k(i)*(khat(1)*x_cc(j) + khat(2)*y_cc(r) + khat(3)*z_cc(l)) + 2._wp*pi*phi
                        velfluc = 2._wp*q*sig*cos(alpha)
                        velfluc = matmul(lmat, velfluc)
                        q_prim_vf(eqn_idx%mom%beg)%sf(j, r, l) = q_prim_vf(eqn_idx%mom%beg)%sf(j, r, l) + velfluc(1)
                        q_prim_vf(eqn_idx%mom%beg + 1)%sf(j, r, l) = q_prim_vf(eqn_idx%mom%beg + 1)%sf(j, r, l) + velfluc(2)
                        q_prim_vf(eqn_idx%mom%beg + 2)%sf(j, r, l) = q_prim_vf(eqn_idx%mom%beg + 2)%sf(j, r, l) + velfluc(3)
                    end do
                end do
            end do
        end do
 
    end subroutine s_perturb_mixlayer
 
    !> Generate deterministic pseudo-random wave vector, polarization, and phase for a perturbation mode.
    subroutine s_generate_random_perturbation(khat, xi, phi, ik, yloc)
 
        integer, intent(in)                 :: ik
        real(wp), intent(in)                :: yloc
        real(wp), dimension(3), intent(out) :: khat, xi
        real(wp), intent(out)               :: phi
        real(wp)                            :: theta, eta
        integer                             :: seed, kfac, yfac
 
        kfac = ik*amplifier
        yfac = nint((sin(yloc) + 1._wp)*amplifier)
        seed = nint(0.5_wp*modmul(kfac) + 0.5_wp*modmul(yfac))
 
        call s_prng(theta, seed)
        call s_prng(eta, seed)
        khat = f_unit_vector(theta, eta)
 
        call s_prng(theta, seed)
        call s_prng(eta, seed)
        xi = f_unit_vector(theta, eta)
 
        call s_prng(phi, seed)
 
    end subroutine s_generate_random_perturbation
 
    !> Generate a unit vector uniformly distributed on the sphere from two random parameters.
    function f_unit_vector(theta, eta) result(vec)
 
        real(wp), intent(in)   :: theta, eta
        real(wp)               :: zeta, xi
        real(wp), dimension(3) :: vec
 
        xi = 2._wp*pi*theta
        zeta = acos(2._wp*eta - 1._wp)
        vec(1) = sin(zeta)*cos(xi)
        vec(2) = sin(zeta)*sin(xi)
        vec(3) = cos(zeta)
 
    end function f_unit_vector
 
    !> Generate a pseudo-random number between 0 and 1 using a linear congruential generator.
    subroutine s_prng(var, seed)
 
        integer, intent(inout) :: seed
        real(wp), intent(out)  :: var
 
        seed = mod(modmul(seed), modulus)
        var = seed/real(modulus, wp)
 
    end subroutine s_prng
 
    !> Compute a modular multiplication step for the linear congruential pseudo-random number generator.
    function modmul(a) result(val)
 
        integer, intent(in) :: a
        integer             :: val
        real(wp)            :: x, y
 
        x = (multiplier/real(modulus, wp))*a + (increment/real(modulus, wp))
        y = nint((x - floor(x))*decimal_trim)/decimal_trim
        val = nint(y*modulus)
 
    end function modmul
 
    !> Deallocate the temporary primitive variable array used by elliptic smoothing.
    impure subroutine s_finalize_perturbation_module()
 
        if (elliptic_smoothing) then
            deallocate (q_prim_temp)
        end if
 
    end subroutine s_finalize_perturbation_module
 
end module m_perturbation