From 328c807bd00b89b19633a309d46628b0f7bc4788 Mon Sep 17 00:00:00 2001 From: Spencer Bryngelson Date: Sat, 14 Jan 2023 12:10:22 -0500 Subject: [PATCH] remove surface tension --- src/post_process/m_derived_variables.f90 | 179 ---- src/simulation/m_data_output.fpp | 909 +----------------- src/simulation/m_derived_variables.f90 | 1106 +--------------------- src/simulation/m_global_parameters.fpp | 10 - src/simulation/m_mpi_proxy.fpp | 4 - src/simulation/m_rhs.fpp | 2 +- src/simulation/m_start_up.fpp | 19 +- src/simulation/m_time_steppers.fpp | 10 +- toolchain/mfc/run/case_dicts.py | 16 +- 9 files changed, 41 insertions(+), 2214 deletions(-) diff --git a/src/post_process/m_derived_variables.f90 b/src/post_process/m_derived_variables.f90 index 16513f4a71..0c9966387c 100644 --- a/src/post_process/m_derived_variables.f90 +++ b/src/post_process/m_derived_variables.f90 @@ -21,14 +21,12 @@ module m_derived_variables private; public :: s_initialize_derived_variables_module, & s_compute_finite_difference_coefficients, & - s_derive_unadvected_volume_fraction, & s_derive_specific_heat_ratio, & s_derive_liquid_stiffness, & s_derive_sound_speed, & s_derive_flux_limiter, & s_derive_vorticity_component, & s_derive_numerical_schlieren_function, & - s_derive_curvature, & s_finalize_derived_variables_module real(kind(0d0)), allocatable, dimension(:, :, :) :: gm_rho_sf !< @@ -168,43 +166,6 @@ subroutine s_compute_finite_difference_coefficients(q, offset_s, & end subroutine s_compute_finite_difference_coefficients ! -------------- - !> This subroutine is used together with the volume fraction - !! model and when called upon, it computes the values of the - !! unadvected volume fraction from the inputted conservative - !! variables, q_cons_vf. The calculated values are stored in - !! the derived flow quantity storage variable, q_sf. - !! @param q_cons_vf Conservative variables - !! @param q_sf Unadvected volume fraction - subroutine s_derive_unadvected_volume_fraction(q_cons_vf, q_sf) ! ------ - - type(scalar_field), & - dimension(sys_size), & - intent(IN) :: q_cons_vf - - real(kind(0d0)), & - dimension(-offset_x%beg:m + offset_x%end, & - -offset_y%beg:n + offset_y%end, & - -offset_z%beg:p + offset_z%end), & - intent(INOUT) :: q_sf - - integer :: i, j, k, l !< Generic loop iterators - - ! Computing unadvected volume fraction from conservative variables - do l = -offset_z%beg, p + offset_z%end - do k = -offset_y%beg, n + offset_y%end - do j = -offset_x%beg, m + offset_x%end - - q_sf(j, k, l) = 1d0 - - do i = adv_idx%beg, adv_idx%end - q_sf(j, k, l) = q_sf(j, k, l) - q_cons_vf(i)%sf(j, k, l) - end do - - end do - end do - end do - - end subroutine s_derive_unadvected_volume_fraction ! ------------------- !> This subroutine receives as input the specific heat ratio !! function, gamma_sf, and derives from it the specific heat @@ -438,146 +399,6 @@ subroutine s_derive_flux_limiter(i, q_prim_vf, q_sf) ! ----------------- end do end subroutine s_derive_flux_limiter ! --------------------------------- - !> Computes the local curvatures - !! @param i Fluid indicator - !! @param q_prim_vf Primitive variables - !! @param q_sf Curvature - subroutine s_derive_curvature(i, q_prim_vf, q_sf) ! -------------------- - - integer, intent(IN) :: i - - type(scalar_field), dimension(sys_size), intent(IN) :: q_prim_vf - - real(kind(0d0)), dimension(-offset_x%beg:m + offset_x%end, & - -offset_y%beg:n + offset_y%end, & - -offset_z%beg:p + offset_z%end), & - intent(INOUT) :: q_sf - - real(kind(0d0)), allocatable, dimension(:, :, :) :: alpha_unadv_sf - - real(kind(0d0)), allocatable, dimension(:, :) :: A - real(kind(0d0)), allocatable, dimension(:) :: sol, b, AA - real(kind(0d0)) :: norm - real(kind(0d0)) :: xloc, yloc, zloc - integer :: ndim - integer :: j, k, l, jj, kk, ll, i1, i2 - integer :: stencil_j_min, stencil_j_max - integer :: stencil_k_min, stencil_k_max - integer :: stencil_l_min, stencil_l_max - type(int_bounds_info) :: iz1 - - if (p > 0) then ! 3D simulation - allocate (A(10, 10)) - allocate (sol(10)); allocate (b(10)); allocate (AA(10)) - else ! 2D simulation - allocate (A(6, 6)) - allocate (sol(6)); allocate (b(6)); allocate (AA(6)) - end if - - if ((i == num_fluids) .and. (adv_alphan .neqv. .true.)) then - allocate (alpha_unadv_sf(-offset_x%beg:m + offset_x%end, & - -offset_y%beg:n + offset_y%end, & - -offset_z%beg:p + offset_z%end)) - call s_derive_unadvected_volume_fraction(q_prim_vf, alpha_unadv_sf) - end if - - if (p > 0) then - iz1%beg = -offset_z%beg; iz1%end = p + offset_z%end - else - iz1%beg = -1; iz1%end = 1 - end if - - ! Parabolic fitting - ndim = size(sol, 1) - - do l = iz1%beg + 1, iz1%end - 1 - do k = -offset_y%beg + 1, n + offset_y%end - 1 - do j = -offset_x%beg + 1, m + offset_x%end - 1 - A(:, :) = 0d0 - b(:) = 0d0 - sol(:) = 0d0 - AA(:) = 0d0 - - stencil_j_min = j - 1; stencil_j_max = j + 1 - stencil_k_min = k - 1; stencil_k_max = k + 1 - if (p > 0) then - stencil_l_min = l - 1; stencil_l_max = l + 1 - else - stencil_l_min = 0; stencil_l_max = 0 - end if - - do ll = stencil_l_min, stencil_l_max - do kk = stencil_k_min, stencil_k_max - do jj = stencil_j_min, stencil_j_max - - ! Ignore corner points in 3D stencil - if ((p > 0) .and. (abs(jj - j) == 1) .and. (abs(kk - k) == 1) .and. (abs(ll - l) == 1)) cycle - - ! Find distance between cell centers - xloc = x_cc(jj) - x_cc(j) - yloc = y_cc(kk) - y_cc(k) - if (p > 0) then - zloc = z_cc(ll) - z_cc(l) - end if - - ! Compute operator - AA(1) = 1d0 - AA(2) = xloc - AA(3) = yloc - AA(4) = 5d-1*xloc**2d0 - AA(5) = 5d-1*yloc**2d0 - AA(6) = xloc*yloc - if (p > 0) then - AA(7) = zloc - AA(8) = 5d-1*zloc**2d0 - AA(9) = yloc*zloc - AA(10) = zloc*xloc - end if - - do i1 = 1, ndim - do i2 = 1, ndim - A(i1, i2) = A(i1, i2) + AA(i1)*AA(i2) - end do - end do - - ! Form RHS vector - do i1 = 1, ndim - if ((i == num_fluids) .and. (adv_alphan .neqv. .true.)) then - b(i1) = b(i1) + alpha_unadv_sf(jj, kk, ll)*AA(i1) - else - b(i1) = b(i1) + q_prim_vf(E_idx + i)%sf(jj, kk, ll)*AA(i1) - end if - end do - end do - end do - end do - - call s_solve_linear_system(A, b, sol, ndim) - - if (p > 0) then - norm = sqrt(sol(2)**2d0 + sol(3)**2d0 + sol(7)**2d0) - else - norm = sqrt(sol(2)**2d0 + sol(3)**2d0) - end if - - if (p > 0) then - q_sf(j, k, l) = -(sol(2)**2d0*sol(5) - 2d0*sol(2)*sol(3)*sol(6) + sol(2)**2d0*sol(8) + & - sol(3)**2d0*sol(8) - 2d0*sol(3)*sol(7)*sol(9) + sol(3)**2d0*sol(4) + & - sol(7)**2d0*sol(4) - 2d0*sol(7)*sol(2)*sol(10) + sol(7)**2d0*sol(5))/ & - max(norm, sgm_eps)**3d0 - else - q_sf(j, k, l) = -(sol(2)**2d0*sol(5) - 2d0*sol(2)*sol(3)*sol(6) + sol(3)**2d0*sol(4))/ & - max(norm, sgm_eps)**3d0 - end if - - end do - end do - end do - - deallocate (A, sol, b, AA) - if ((i == num_fluids) .and. (adv_alphan .neqv. .true.)) deallocate (alpha_unadv_sf) - - end subroutine s_derive_curvature ! ------------------------------------ !> Computes the solution to the linear system Ax=b w/ sol = x !! @param A Input matrix diff --git a/src/simulation/m_data_output.fpp b/src/simulation/m_data_output.fpp index e7fcfc5660..18dfa8d549 100644 --- a/src/simulation/m_data_output.fpp +++ b/src/simulation/m_data_output.fpp @@ -29,19 +29,13 @@ module m_data_output private; public :: s_initialize_data_output_module, & s_open_run_time_information_file, & - s_open_com_files, & - s_open_cb_files, & s_open_probe_files, & s_write_run_time_information, & s_write_data_files, & s_write_serial_data_files, & s_write_parallel_data_files, & - s_write_com_files, & - s_write_cb_files, & s_write_probe_files, & s_close_run_time_information_file, & - s_close_com_files, & - s_close_cb_files, & s_close_probe_files, & s_finalize_data_output_module @@ -85,19 +79,11 @@ module m_data_output real(kind(0d0)) :: Rc_min !< Rc criterion maximum !> @} - ! @name Generic storage for flow variable(s) that are to be written to CoM data file + ! @name Variables for computing acceleration !> @{ real(kind(0d0)), public, allocatable, dimension(:, :, :) :: accel_mag - real(kind(0d0)), public, allocatable, dimension(:, :) :: q_com - real(kind(0d0)), public, allocatable, dimension(:, :, :) :: moments - real(kind(0d0)), public, allocatable, dimension(:, :) :: cb_mass - real(kind(0d0)), public, allocatable, dimension(:, :, :) :: bounds - real(kind(0d0)), public, allocatable, dimension(:, :) :: cntrline real(kind(0d0)), public, allocatable, dimension(:, :, :) :: x_accel, y_accel, z_accel - type(scalar_field), allocatable, dimension(:) :: grad_x_vf, grad_y_vf, grad_z_vf, norm_vf - real(kind(0d0)), target, allocatable, dimension(:, :, :) :: energy - - + !> @} procedure(s_write_abstract_data_files), pointer :: s_write_data_files => null() @@ -169,131 +155,6 @@ contains end subroutine s_open_run_time_information_file ! ---------------------- - !> This opens a formatted data file where the root processor - !! can write out the CoM information - subroutine s_open_com_files() ! ---------------------------------------- - - character(LEN=path_len + 3*name_len) :: file_path !< - !! Relative path to the CoM file in the case directory - - integer :: i !< Generic loop iterator - - do i = 1, num_fluids - if (com_wrt(i)) then - - ! Generating the relative path to the CoM data file - write (file_path, '(A,I0,A)') '/fluid', i, '_com.dat' - file_path = trim(case_dir)//trim(file_path) - - ! Creating the formatted data file and setting up its - ! structure - open (i + 10, FILE=trim(file_path), & - FORM='formatted', & - POSITION='append', & - STATUS='unknown') - - if (n == 0) then - if (any(moment_order /= dflt_int)) then - write (i + 10, '(A)') '=== Non-Dimensional Time '// & - '=== Total Mass '// & - '=== x-loc '// & - '=== x-vel '// & - '=== x-accel '// & - '=== Higher Moments ===' - else - write (i + 10, '(A)') '=== Non-Dimensional Time '// & - '=== Total Mass '// & - '=== x-loc '// & - '=== x-vel '// & - '=== x-accel ===' - end if - elseif (p == 0) then - if (any(moment_order /= dflt_int)) then - write (i + 10, '(A)') '=== Non-Dimensional Time '// & - '=== Total Mass '// & - '=== x-loc '// & - '=== y-loc '// & - '=== x-vel '// & - '=== y-vel '// & - '=== x-accel '// & - '=== y-accel '// & - '=== Higher Moments ===' - else - write (i + 10, '(A)') '=== Non-Dimensional Time '// & - '=== Total Mass '// & - '=== x-loc '// & - '=== y-loc '// & - '=== x-vel '// & - '=== y-vel '// & - '=== x-accel '// & - '=== y-accel ===' - end if - else - if (any(moment_order /= dflt_int)) then - write (i + 10, '(A)') '=== Non-Dimensional Time '// & - '=== Total Mass '// & - '=== x-loc '// & - '=== y-loc '// & - '=== z-loc '// & - '=== x-vel '// & - '=== y-vel '// & - '=== z-vel '// & - '=== x-accel '// & - '=== y-accel '// & - '=== z-accel '// & - '=== Higher Moments ===' - else - write (i + 10, '(A)') '=== Non-Dimensional Time '// & - '=== Total Mass '// & - '=== x-loc '// & - '=== y-loc '// & - '=== z-loc '// & - '=== x-vel '// & - '=== y-vel '// & - '=== z-vel '// & - '=== x-accel '// & - '=== y-accel '// & - '=== z-accel ===' - end if - end if - end if - end do - - end subroutine s_open_com_files ! -------------------------------------- - - !> This opens a formatted data file where the root processor - !! can write out the coherent body (cb) information - subroutine s_open_cb_files() ! ---------------------------------------- - - character(LEN=path_len + 3*name_len) :: file_path !< - !! Relative path to the cb file in the case directory - - integer :: i !< Generic loop iterator - - do i = 1, num_fluids - if (cb_wrt(i)) then - - ! Generating the relative path to the cb data file - write (file_path, '(A,I0,A)') '/fluid', i, '_cb.dat' - file_path = trim(case_dir)//trim(file_path) - - ! Creating the formatted data file and setting up its - ! structure - open (i + 20, FILE=trim(file_path), & - FORM='formatted', & - POSITION='append', & - STATUS='unknown') - - write (i + 20, '(A)') '=== Non-Dimensional Time '// & - '=== Coherent Body Mass '// & - '=== Coherent Body Area '// & - '=== Bounds '// & - '=== Centerline ===' - end if - end do - - end subroutine s_open_cb_files ! -------------------------------------- - !> This opens a formatted data file where the root processor !! can write out flow probe information subroutine s_open_probe_files() ! -------------------------------------- @@ -914,718 +775,6 @@ contains end subroutine s_write_parallel_data_files ! --------------------------- - !> This writes a formatted data file where the root processor - !! can write out the CoM information - !! @param t_step Current time-step - !! @param q_com Center of mass information - !! @param moments Higher moment information - subroutine s_write_com_files(t_step, q_com, moments) ! ------------------- - - integer, intent(IN) :: t_step - real(kind(0d0)), dimension(num_fluids, 10), intent(IN) :: q_com - real(kind(0d0)), dimension(num_fluids, 2, 5), intent(IN) :: moments - - integer :: i !< Generic loop iterator - real(kind(0d0)) :: nondim_time !< Non-dimensional time - - ! Non-dimensional time calculation - if (t_step_old /= dflt_int) then - nondim_time = real(t_step + t_step_old, kind(0d0))*dt - else - nondim_time = real(t_step, kind(0d0))*dt - end if - - if (n == 0) then ! 1D simulation - do i = 1, num_fluids ! Loop through fluids - if (com_wrt(i)) then ! Writing out CoM data - if (proc_rank == 0) then - write (i + 10, '(6X,F12.6,F24.8,F24.8,F24.8,F24.8)') & - nondim_time, & - q_com(i, 1), & - q_com(i, 2), & - q_com(i, 5), & - q_com(i, 8) - - end if - end if - end do - elseif (p == 0) then ! 2D simulation - do i = 1, num_fluids ! Loop through fluids - if (com_wrt(i)) then ! Writing out CoM data - if (proc_rank == 0) then - if (moment_order(1) == dflt_int) then - write (i + 10, '(6X,F12.6,F24.8,F24.8,F24.8,F24.8,'// & - 'F24.8,F24.8,F24.8)') & - nondim_time, & - q_com(i, 1), & - q_com(i, 2), & - q_com(i, 3), & - q_com(i, 5), & - q_com(i, 6), & - q_com(i, 8), & - q_com(i, 9) - elseif (moment_order(2) == dflt_int) then - write (i + 10, '(6X,F12.6,F24.8,F24.8,F24.8,F24.8,'// & - 'F24.8,F24.8,F24.8,E24.8)') & - nondim_time, & - q_com(i, 1), & - q_com(i, 2), & - q_com(i, 3), & - q_com(i, 5), & - q_com(i, 6), & - q_com(i, 8), & - q_com(i, 9), & - moments(i, 1, 1) - elseif (moment_order(3) == dflt_int) then - write (i + 10, '(6X,F12.6,F24.8,F24.8,F24.8,F24.8,'// & - 'F24.8,F24.8,F24.8,E24.8,'// & - 'E24.8)') & - nondim_time, & - q_com(i, 1), & - q_com(i, 2), & - q_com(i, 3), & - q_com(i, 5), & - q_com(i, 6), & - q_com(i, 8), & - q_com(i, 9), & - moments(i, 1, 1), & - moments(i, 1, 2) - elseif (moment_order(4) == dflt_int) then - write (i + 10, '(6X,F12.6,F24.8,F24.8,F24.8,F24.8,'// & - 'F24.8,F24.8,F24.8,E24.8,'// & - 'E24.8,E24.8)') & - nondim_time, & - q_com(i, 1), & - q_com(i, 2), & - q_com(i, 3), & - q_com(i, 5), & - q_com(i, 6), & - q_com(i, 8), & - q_com(i, 9), & - moments(i, 1, 1), & - moments(i, 1, 2), & - moments(i, 1, 3) - elseif (moment_order(5) == dflt_int) then - write (i + 10, '(6X,F12.6,F24.8,F24.8,F24.8,F24.8,'// & - 'F24.8,F24.8,F24.8,E24.8,'// & - 'E24.8,E24.8,E24.8)') & - nondim_time, & - q_com(i, 1), & - q_com(i, 2), & - q_com(i, 3), & - q_com(i, 5), & - q_com(i, 6), & - q_com(i, 8), & - q_com(i, 9), & - moments(i, 1, 1), & - moments(i, 1, 2), & - moments(i, 1, 3), & - moments(i, 1, 4) - else - write (i + 10, '(6X,F12.6,F24.8,F24.8,F24.8,F24.8,'// & - 'F24.8,F24.8,F24.8,E24.8,'// & - 'E24.8,E24.8,E24.8,E24.8)') & - nondim_time, & - q_com(i, 1), & - q_com(i, 2), & - q_com(i, 3), & - q_com(i, 5), & - q_com(i, 6), & - q_com(i, 8), & - q_com(i, 9), & - moments(i, 1, 1), & - moments(i, 1, 2), & - moments(i, 1, 3), & - moments(i, 1, 4), & - moments(i, 1, 5) - end if - end if - end if - end do - else ! 3D simulation - do i = 1, num_fluids ! Loop through fluids - if (com_wrt(i)) then ! Writing out CoM data - if (proc_rank == 0) then - if (moment_order(1) == dflt_int) then - write (i + 10, '(6X,F12.6,F24.8,F24.8,F24.8,F24.8,'// & - 'F24.8,F24.8,F24.8,F24.8,'// & - 'F24.8,F24.8)') & - nondim_time, & - q_com(i, 1), & - q_com(i, 2), & - q_com(i, 3), & - q_com(i, 4), & - q_com(i, 5), & - q_com(i, 6), & - q_com(i, 7), & - q_com(i, 8), & - q_com(i, 9), & - q_com(i, 10) - elseif (moment_order(2) == dflt_int) then - write (i + 10, '(6X,F12.6,F24.8,F24.8,F24.8,F24.8,'// & - 'F24.8,F24.8,F24.8,F24.8,'// & - 'F24.8,F24.8,F24.8,F24.8)') & - nondim_time, & - q_com(i, 1), & - q_com(i, 2), & - q_com(i, 3), & - q_com(i, 4), & - q_com(i, 5), & - q_com(i, 6), & - q_com(i, 7), & - q_com(i, 8), & - q_com(i, 9), & - q_com(i, 10), & - moments(i, 1, 1), & - moments(i, 2, 1) - elseif (moment_order(3) == dflt_int) then - write (i + 10, '(6X,F12.6,F24.8,F24.8,F24.8,F24.8,'// & - 'F24.8,F24.8,F24.8,F24.8,'// & - 'F24.8,F24.8,F24.8,F24.8,'// & - 'F24.8,F24.8)') & - nondim_time, & - q_com(i, 1), & - q_com(i, 2), & - q_com(i, 3), & - q_com(i, 4), & - q_com(i, 5), & - q_com(i, 6), & - q_com(i, 7), & - q_com(i, 8), & - q_com(i, 9), & - q_com(i, 10), & - moments(i, 1, 1), & - moments(i, 1, 2), & - moments(i, 2, 1), & - moments(i, 2, 2) - elseif (moment_order(4) == dflt_int) then - write (i + 10, '(6X,F12.6,F24.8,F24.8,F24.8,F24.8,'// & - 'F24.8,F24.8,F24.8,F24.8,'// & - 'F24.8,F24.8,F24.8,F24.8,'// & - 'F24.8,F24.8,F24.8,F24.8)') & - nondim_time, & - q_com(i, 1), & - q_com(i, 2), & - q_com(i, 3), & - q_com(i, 4), & - q_com(i, 5), & - q_com(i, 6), & - q_com(i, 7), & - q_com(i, 8), & - q_com(i, 9), & - q_com(i, 10), & - moments(i, 1, 1), & - moments(i, 1, 2), & - moments(i, 1, 3), & - moments(i, 2, 1), & - moments(i, 2, 2), & - moments(i, 2, 3) - elseif (moment_order(5) == dflt_int) then - write (i + 10, '(6X,F12.6,F24.8,F24.8,F24.8,F24.8,'// & - 'F24.8,F24.8,F24.8,F24.8,'// & - 'F24.8,F24.8,F24.8,F24.8,'// & - 'F24.8,F24.8,F24.8,F24.8,'// & - 'F24.8,F24.8)') & - nondim_time, & - q_com(i, 1), & - q_com(i, 2), & - q_com(i, 3), & - q_com(i, 4), & - q_com(i, 5), & - q_com(i, 6), & - q_com(i, 7), & - q_com(i, 8), & - q_com(i, 9), & - q_com(i, 10), & - moments(i, 1, 1), & - moments(i, 1, 2), & - moments(i, 1, 3), & - moments(i, 1, 4), & - moments(i, 2, 1), & - moments(i, 2, 2), & - moments(i, 2, 3), & - moments(i, 2, 4) - else - write (i + 10, '(6X,F12.6,F24.8,F24.8,F24.8,F24.8,'// & - 'F24.8,F24.8,F24.8,F24.8,'// & - 'F24.8,F24.8,F24.8,F24.8,'// & - 'F24.8,F24.8,F24.8,F24.8,'// & - 'F24.8,F24.8,F24.8,F24.8)') & - nondim_time, & - q_com(i, 1), & - q_com(i, 2), & - q_com(i, 3), & - q_com(i, 4), & - q_com(i, 5), & - q_com(i, 6), & - q_com(i, 7), & - q_com(i, 8), & - q_com(i, 9), & - q_com(i, 10), & - moments(i, 1, 1), & - moments(i, 1, 2), & - moments(i, 1, 3), & - moments(i, 1, 4), & - moments(i, 1, 5), & - moments(i, 2, 1), & - moments(i, 2, 2), & - moments(i, 2, 3), & - moments(i, 2, 4), & - moments(i, 2, 5) - end if - end if - end if - end do - end if - - end subroutine s_write_com_files ! ------------------------------------- - - !> The goal of this subroutine is to output coherent body information. - !! @param t_step Current time-step - !! @param cb_mass Coherent body mass - !! @param bounds Coherent body boundary - !! @param cntrline Coherent body center line - subroutine s_write_cb_files(t_step, cb_mass, bounds, cntrline) ! ---------- - - integer, intent(IN) :: t_step - real(kind(0d0)), dimension(num_fluids, 10), intent(IN) :: cb_mass - real(kind(0d0)), dimension(num_fluids, 5, 6), intent(IN) :: bounds - real(kind(0d0)), dimension(num_fluids, 5), intent(IN) :: cntrline - - integer :: i !< Generic loop iterator - real(kind(0d0)) :: nondim_time !< Non-dimensional time - - ! Non-dimensional time calculation - if (t_step_old /= dflt_int) then - nondim_time = real(t_step + t_step_old, kind(0d0))*dt - else - nondim_time = real(t_step, kind(0d0))*dt - end if - - do i = 1, num_fluids ! Loop through fluids - if (cb_wrt(i)) then ! Writing out CoM data - if (proc_rank == 0) then - if (n == 0) then ! 1D simulation - if (threshold_mf(2) == dflt_real) then - write (i + 20, '(6X,F12.6,F24.8,F24.8,F24.8,F24.8)') & - nondim_time, & - cb_mass(i, 1), & - cb_mass(i, 6), & - bounds(i, 1, 1), & - bounds(i, 1, 2) - elseif (threshold_mf(3) == dflt_real) then - write (i + 20, '(6X,F12.6,F24.8,F24.8,F24.8,F24.8,'// & - 'F24.8,F24.8,F24.8,F24.8)') & - nondim_time, & - cb_mass(i, 1), & - cb_mass(i, 2), & - cb_mass(i, 6), & - cb_mass(i, 7), & - bounds(i, 1, 1), & - bounds(i, 2, 1), & - bounds(i, 1, 2), & - bounds(i, 2, 2) - elseif (threshold_mf(4) == dflt_real) then - write (i + 20, '(6X,F12.6,F24.8,F24.8,F24.8,F24.8,'// & - 'F24.8,F24.8,F24.8,F24.8,'// & - 'F24.8,F24.8,F24.8,F24.8)') & - nondim_time, & - cb_mass(i, 1), & - cb_mass(i, 2), & - cb_mass(i, 3), & - cb_mass(i, 6), & - cb_mass(i, 7), & - cb_mass(i, 8), & - bounds(i, 1, 1), & - bounds(i, 2, 1), & - bounds(i, 3, 1), & - bounds(i, 1, 2), & - bounds(i, 2, 2), & - bounds(i, 3, 2) - elseif (threshold_mf(5) == dflt_real) then - write (i + 20, '(6X,F12.6,F24.8,F24.8,F24.8,F24.8,'// & - 'F24.8,F24.8,F24.8,F24.8,'// & - 'F24.8,F24.8,F24.8,F24.8,'// & - 'F24.8,F24.8,F24.8,F24.8)') & - nondim_time, & - cb_mass(i, 1), & - cb_mass(i, 2), & - cb_mass(i, 3), & - cb_mass(i, 4), & - cb_mass(i, 6), & - cb_mass(i, 7), & - cb_mass(i, 8), & - cb_mass(i, 9), & - bounds(i, 1, 1), & - bounds(i, 2, 1), & - bounds(i, 3, 1), & - bounds(i, 4, 1), & - bounds(i, 1, 2), & - bounds(i, 2, 2), & - bounds(i, 3, 2), & - bounds(i, 4, 2) - else - write (i + 20, '(6X,F12.6,F24.8,F24.8,F24.8,F24.8,'// & - 'F24.8,F24.8,F24.8,F24.8,'// & - 'F24.8,F24.8,F24.8,F24.8,'// & - 'F24.8,F24.8,F24.8,F24.8,'// & - 'F24.8,F24.8,F24.8,F24.8)') & - nondim_time, & - cb_mass(i, 1), & - cb_mass(i, 2), & - cb_mass(i, 3), & - cb_mass(i, 4), & - cb_mass(i, 5), & - cb_mass(i, 6), & - cb_mass(i, 7), & - cb_mass(i, 8), & - cb_mass(i, 9), & - cb_mass(i, 10), & - bounds(i, 1, 1), & - bounds(i, 2, 1), & - bounds(i, 3, 1), & - bounds(i, 4, 1), & - bounds(i, 5, 1), & - bounds(i, 1, 2), & - bounds(i, 2, 2), & - bounds(i, 3, 2), & - bounds(i, 4, 2), & - bounds(i, 5, 2) - end if - elseif (p == 0) then ! 2D simulation - if (threshold_mf(2) == dflt_real) then - write (i + 20, '(6X,F12.6,F24.8,F24.8,F24.8,F24.8,'// & - 'F24.8,F24.8,F24.8)') & - nondim_time, & - cb_mass(i, 1), & - cb_mass(i, 6), & - bounds(i, 1, 1), & - bounds(i, 1, 2), & - bounds(i, 1, 3), & - bounds(i, 1, 4), & - cntrline(i, 1) - elseif (threshold_mf(3) == dflt_real) then - write (i + 20, '(6X,F12.6,F24.8,F24.8,F24.8,F24.8,'// & - 'F24.8,F24.8,F24.8,F24.8,'// & - 'F24.8,F24.8,F24.8,F24.8,'// & - 'F24.8,F24.8)') & - nondim_time, & - cb_mass(i, 1), & - cb_mass(i, 2), & - cb_mass(i, 6), & - cb_mass(i, 7), & - bounds(i, 1, 1), & - bounds(i, 2, 1), & - bounds(i, 1, 2), & - bounds(i, 2, 2), & - bounds(i, 1, 3), & - bounds(i, 2, 3), & - bounds(i, 1, 4), & - bounds(i, 2, 4), & - cntrline(i, 1), & - cntrline(i, 2) - elseif (threshold_mf(4) == dflt_real) then - write (i + 20, '(6X,F12.6,F24.8,F24.8,F24.8,F24.8,'// & - 'F24.8,F24.8,F24.8,F24.8,'// & - 'F24.8,F24.8,F24.8,F24.8,'// & - 'F24.8,F24.8,F24.8,F24.8,'// & - 'F24.8,F24.8,F24.8,F24.8,'// & - 'F24.8)') & - nondim_time, & - cb_mass(i, 1), & - cb_mass(i, 2), & - cb_mass(i, 3), & - cb_mass(i, 6), & - cb_mass(i, 7), & - cb_mass(i, 8), & - bounds(i, 1, 1), & - bounds(i, 2, 1), & - bounds(i, 3, 1), & - bounds(i, 1, 2), & - bounds(i, 2, 2), & - bounds(i, 3, 2), & - bounds(i, 1, 3), & - bounds(i, 2, 3), & - bounds(i, 3, 3), & - bounds(i, 1, 4), & - bounds(i, 2, 4), & - bounds(i, 3, 4), & - cntrline(i, 1), & - cntrline(i, 2), & - cntrline(i, 3) - elseif (threshold_mf(5) == dflt_real) then - write (i + 20, '(6X,F12.6,F24.8,F24.8,F24.8,F24.8,'// & - 'F24.8,F24.8,F24.8,F24.8,'// & - 'F24.8,F24.8,F24.8,F24.8,'// & - 'F24.8,F24.8,F24.8,F24.8,'// & - 'F24.8,F24.8,F24.8,F24.8,'// & - 'F24.8,F24.8,F24.8,F24.8,'// & - 'F24.8,F24.8,F24.8,F24.8)') & - nondim_time, & - cb_mass(i, 1), & - cb_mass(i, 2), & - cb_mass(i, 3), & - cb_mass(i, 4), & - cb_mass(i, 6), & - cb_mass(i, 7), & - cb_mass(i, 8), & - cb_mass(i, 9), & - bounds(i, 1, 1), & - bounds(i, 2, 1), & - bounds(i, 3, 1), & - bounds(i, 4, 1), & - bounds(i, 1, 2), & - bounds(i, 2, 2), & - bounds(i, 3, 2), & - bounds(i, 4, 2), & - bounds(i, 1, 3), & - bounds(i, 2, 3), & - bounds(i, 3, 3), & - bounds(i, 4, 3), & - bounds(i, 1, 4), & - bounds(i, 2, 4), & - bounds(i, 3, 4), & - bounds(i, 4, 4), & - cntrline(i, 1), & - cntrline(i, 2), & - cntrline(i, 3), & - cntrline(i, 4) - else - write (i + 20, '(6X,F12.6,F24.8,F24.8,F24.8,F24.8,'// & - 'F24.8,F24.8,F24.8,F24.8,'// & - 'F24.8,F24.8,F24.8,F24.8,'// & - 'F24.8,F24.8,F24.8,F24.8,'// & - 'F24.8,F24.8,F24.8,F24.8,'// & - 'F24.8,F24.8,F24.8,F24.8,'// & - 'F24.8,F24.8,F24.8,F24.8,'// & - 'F24.8,F24.8,F24.8,F24.8,'// & - 'F24.8,F24.8,F24.8)') & - nondim_time, & - cb_mass(i, 1), & - cb_mass(i, 2), & - cb_mass(i, 3), & - cb_mass(i, 4), & - cb_mass(i, 5), & - cb_mass(i, 6), & - cb_mass(i, 7), & - cb_mass(i, 8), & - cb_mass(i, 9), & - cb_mass(i, 10), & - bounds(i, 1, 1), & - bounds(i, 2, 1), & - bounds(i, 3, 1), & - bounds(i, 4, 1), & - bounds(i, 5, 1), & - bounds(i, 1, 2), & - bounds(i, 2, 2), & - bounds(i, 3, 2), & - bounds(i, 4, 2), & - bounds(i, 5, 2), & - bounds(i, 1, 3), & - bounds(i, 2, 3), & - bounds(i, 3, 3), & - bounds(i, 4, 3), & - bounds(i, 5, 3), & - bounds(i, 1, 4), & - bounds(i, 2, 4), & - bounds(i, 3, 4), & - bounds(i, 4, 4), & - bounds(i, 5, 4), & - cntrline(i, 1), & - cntrline(i, 2), & - cntrline(i, 3), & - cntrline(i, 4), & - cntrline(i, 5) - end if - else ! 3D simulation - if (threshold_mf(2) == dflt_real) then - write (i + 20, '(6X,F12.6,F24.8,F24.8,F24.8,F24.8,'// & - 'F24.8,F24.8,F24.8,F24.8,'// & - 'F24.8)') & - nondim_time, & - cb_mass(i, 1), & - cb_mass(i, 6), & - bounds(i, 1, 1), & - bounds(i, 1, 2), & - bounds(i, 1, 3), & - bounds(i, 1, 4), & - bounds(i, 1, 5), & - bounds(i, 1, 6), & - cntrline(i, 1) - elseif (threshold_mf(3) == dflt_real) then - write (i + 20, '(6X,F12.6,F24.8,F24.8,F24.8,F24.8,'// & - 'F24.8,F24.8,F24.8,F24.8,'// & - 'F24.8,F24.8,F24.8,F24.8,'// & - 'F24.8,F24.8,F24.8,F24.8,'// & - 'F24.8,F24.8)') & - nondim_time, & - cb_mass(i, 1), & - cb_mass(i, 2), & - cb_mass(i, 6), & - cb_mass(i, 7), & - bounds(i, 1, 1), & - bounds(i, 2, 1), & - bounds(i, 1, 2), & - bounds(i, 2, 2), & - bounds(i, 1, 3), & - bounds(i, 2, 3), & - bounds(i, 1, 4), & - bounds(i, 2, 4), & - bounds(i, 1, 5), & - bounds(i, 2, 5), & - bounds(i, 1, 6), & - bounds(i, 2, 6), & - cntrline(i, 1), & - cntrline(i, 2) - elseif (threshold_mf(4) == dflt_real) then - write (i + 20, '(6X,F12.6,F24.8,F24.8,F24.8,F24.8,'// & - 'F24.8,F24.8,F24.8,F24.8,'// & - 'F24.8,F24.8,F24.8,F24.8,'// & - 'F24.8,F24.8,F24.8,F24.8,'// & - 'F24.8,F24.8,F24.8,F24.8,'// & - 'F24.8,F24.8,F24.8,F24.8,'// & - 'F24.8,F24.8,F24.8)') & - nondim_time, & - cb_mass(i, 1), & - cb_mass(i, 2), & - cb_mass(i, 3), & - cb_mass(i, 6), & - cb_mass(i, 7), & - cb_mass(i, 8), & - bounds(i, 1, 1), & - bounds(i, 2, 1), & - bounds(i, 3, 1), & - bounds(i, 1, 2), & - bounds(i, 2, 2), & - bounds(i, 3, 2), & - bounds(i, 1, 3), & - bounds(i, 2, 3), & - bounds(i, 3, 3), & - bounds(i, 1, 4), & - bounds(i, 2, 4), & - bounds(i, 3, 4), & - bounds(i, 1, 5), & - bounds(i, 2, 5), & - bounds(i, 3, 5), & - bounds(i, 1, 6), & - bounds(i, 2, 6), & - bounds(i, 3, 6), & - cntrline(i, 1), & - cntrline(i, 2), & - cntrline(i, 3) - elseif (threshold_mf(5) == dflt_real) then - write (i + 20, '(6X,F12.6,F24.8,F24.8,F24.8,F24.8,'// & - 'F24.8,F24.8,F24.8,F24.8,'// & - 'F24.8,F24.8,F24.8,F24.8,'// & - 'F24.8,F24.8,F24.8,F24.8,'// & - 'F24.8,F24.8,F24.8,F24.8,'// & - 'F24.8,F24.8,F24.8,F24.8,'// & - 'F24.8,F24.8,F24.8,F24.8,'// & - 'F24.8,F24.8,F24.8,F24.8,'// & - 'F24.8,F24.8,F24.8,F24.8)') & - nondim_time, & - cb_mass(i, 1), & - cb_mass(i, 2), & - cb_mass(i, 3), & - cb_mass(i, 4), & - cb_mass(i, 6), & - cb_mass(i, 7), & - cb_mass(i, 8), & - cb_mass(i, 9), & - bounds(i, 1, 1), & - bounds(i, 2, 1), & - bounds(i, 3, 1), & - bounds(i, 4, 1), & - bounds(i, 1, 2), & - bounds(i, 2, 2), & - bounds(i, 3, 2), & - bounds(i, 4, 2), & - bounds(i, 1, 3), & - bounds(i, 2, 3), & - bounds(i, 3, 3), & - bounds(i, 4, 3), & - bounds(i, 1, 4), & - bounds(i, 2, 4), & - bounds(i, 3, 4), & - bounds(i, 4, 4), & - bounds(i, 1, 5), & - bounds(i, 2, 5), & - bounds(i, 3, 5), & - bounds(i, 4, 5), & - bounds(i, 1, 6), & - bounds(i, 2, 6), & - bounds(i, 3, 6), & - bounds(i, 4, 6), & - cntrline(i, 1), & - cntrline(i, 2), & - cntrline(i, 3), & - cntrline(i, 4) - else - write (i + 20, '(6X,F12.6,F24.8,F24.8,F24.8,F24.8,'// & - 'F24.8,F24.8,F24.8,F24.8,'// & - 'F24.8,F24.8,F24.8,F24.8,'// & - 'F24.8,F24.8,F24.8,F24.8,'// & - 'F24.8,F24.8,F24.8,F24.8,'// & - 'F24.8,F24.8,F24.8,F24.8,'// & - 'F24.8,F24.8,F24.8,F24.8,'// & - 'F24.8,F24.8,F24.8,F24.8,'// & - 'F24.8,F24.8,F24.8,F24.8,'// & - 'F24.8,F24.8,F24.8,F24.8,'// & - 'F24.8,F24.8,F24.8,F24.8,'// & - 'F24.8)') & - nondim_time, & - cb_mass(i, 1), & - cb_mass(i, 2), & - cb_mass(i, 3), & - cb_mass(i, 4), & - cb_mass(i, 5), & - cb_mass(i, 6), & - cb_mass(i, 7), & - cb_mass(i, 8), & - cb_mass(i, 9), & - cb_mass(i, 10), & - bounds(i, 1, 1), & - bounds(i, 2, 1), & - bounds(i, 3, 1), & - bounds(i, 4, 1), & - bounds(i, 5, 1), & - bounds(i, 1, 2), & - bounds(i, 2, 2), & - bounds(i, 3, 2), & - bounds(i, 4, 2), & - bounds(i, 5, 2), & - bounds(i, 1, 3), & - bounds(i, 2, 3), & - bounds(i, 3, 3), & - bounds(i, 4, 3), & - bounds(i, 5, 3), & - bounds(i, 1, 4), & - bounds(i, 2, 4), & - bounds(i, 3, 4), & - bounds(i, 4, 4), & - bounds(i, 5, 4), & - bounds(i, 1, 5), & - bounds(i, 2, 5), & - bounds(i, 3, 5), & - bounds(i, 4, 5), & - bounds(i, 5, 5), & - bounds(i, 1, 6), & - bounds(i, 2, 6), & - bounds(i, 3, 6), & - bounds(i, 4, 6), & - bounds(i, 5, 6), & - cntrline(i, 1), & - cntrline(i, 2), & - cntrline(i, 3), & - cntrline(i, 4), & - cntrline(i, 5) - end if - end if - end if - end if - end do - - end subroutine s_write_cb_files ! ------------------------------------- !> This writes a formatted data file for the flow probe information !! @param t_step Current time-step @@ -2357,32 +1506,6 @@ contains end subroutine s_close_run_time_information_file ! --------------------- - !> Closes communication files - subroutine s_close_com_files() ! --------------------------------------- - - integer :: i !< Generic loop iterator - - do i = 1, num_fluids - if (com_wrt(i)) then - close (i + 10) - end if - end do - - end subroutine s_close_com_files ! ------------------------------------- - - !> Closes coherent body files - subroutine s_close_cb_files() ! --------------------------------------- - - integer :: i !< Generic loop iterator - - do i = 1, num_fluids - if (cb_wrt(i)) then - close (i + 20) - end if - end do - - end subroutine s_close_cb_files ! ------------------------------------- - !> Closes probe files subroutine s_close_probe_files() ! ------------------------------------- @@ -2430,24 +1553,6 @@ contains s_convert_species_to_mixture_variables end if - ! Allocating the generic storage for the flow variable(s) that are - ! going to be written to the CoM data files - if (any(com_wrt)) then - ! num_fluids, mass, x-loc, y-loc, z-loc, x-vel, y-vel, z-vel, x-acc, y-acc, z-acc - allocate (q_com(num_fluids, 10)) - ! num_fluids, 2 lateral directions, 5 higher moment orders - allocate (moments(num_fluids, 2, 5)) - end if - - if (any(cb_wrt)) then - ! num_fluids, mass for 5 threshold mass fractions, area for 5 threshold mass fractions - allocate (cb_mass(num_fluids, 10)) - ! num_fluids, 5 threshold mass fractions, xmin, xmax, ymin, ymax, zmin, zmax - allocate (bounds(num_fluids, 5, 6)) - ! num_fluids, centerline dimension for 5 threshold mass fractions - allocate (cntrline(num_fluids, 5)) - end if - if (probe_wrt) then allocate (accel_mag(0:m, 0:n, 0:p)) allocate (x_accel(0:m, 0:n, 0:p)) @@ -2478,16 +1583,6 @@ contains @:DEALLOCATE(vcfl_sf, Rc_sf) end if - ! Deallocating the storage employed for the flow variables that - ! were written to CoM and probe files - if (any(com_wrt)) then - deallocate(q_com, moments) - end if - - if (any(cb_wrt)) then - deallocate(cb_mass, bounds, cntrline) - end if - if (probe_wrt) then deallocate(accel_mag, x_accel) if (n > 0) then diff --git a/src/simulation/m_derived_variables.f90 b/src/simulation/m_derived_variables.f90 index 08b91513d6..253eb04ec3 100644 --- a/src/simulation/m_derived_variables.f90 +++ b/src/simulation/m_derived_variables.f90 @@ -71,12 +71,6 @@ subroutine s_initialize_derived_variables() ! ----------------------------- ! Opening and writing header of CoM and flow probe files if (proc_rank == 0) then - if (any(com_wrt)) then - call s_open_com_files() - end if - if (any(cb_wrt)) then - call s_open_cb_files() - end if if (probe_wrt) then call s_open_probe_files() end if @@ -103,65 +97,46 @@ subroutine s_compute_derived_variables(t_step) ! ----------------------- integer :: i, j, k !< Generic loop iterators - ! IF ((ANY(com_wrt) .OR. ANY(cb_wrt) .OR. probe_wrt) .AND. (t_step > t_step_start + 2)) THEN - if ((any(com_wrt) .or. any(cb_wrt) .or. probe_wrt)) then - if (any(com_wrt)) then - call s_derive_center_of_mass(q_prim_ts(0)%vf, & - q_prim_ts(1)%vf, & - q_prim_ts(2)%vf, & - q_prim_ts(3)%vf, & - q_com) - call s_derive_higher_moments(q_prim_ts(0)%vf, moments) - call s_write_com_files(t_step, q_com, moments) - end if - - if (any(cb_wrt)) then - call s_derive_fluid_bounds(q_prim_ts(0)%vf, bounds) - call s_derive_coherent_body(q_prim_ts(0)%vf, cb_mass) - call s_derive_centerline(q_prim_ts(0)%vf, cntrline) - call s_write_cb_files(t_step, cb_mass, bounds, cntrline) - end if - - if (probe_wrt) then - call s_derive_acceleration_component(1, q_prim_ts(0)%vf, & + ! IF ( probe_wrt .AND. (t_step > t_step_start + 2)) THEN + if (probe_wrt) then + call s_derive_acceleration_component(1, q_prim_ts(0)%vf, & + q_prim_ts(1)%vf, & + q_prim_ts(2)%vf, & + q_prim_ts(3)%vf, & + x_accel) + if (n > 0) then + call s_derive_acceleration_component(2, q_prim_ts(0)%vf, & q_prim_ts(1)%vf, & q_prim_ts(2)%vf, & q_prim_ts(3)%vf, & - x_accel) - if (n > 0) then - call s_derive_acceleration_component(2, q_prim_ts(0)%vf, & + y_accel) + if (p > 0) then + call s_derive_acceleration_component(3, q_prim_ts(0)%vf, & q_prim_ts(1)%vf, & q_prim_ts(2)%vf, & q_prim_ts(3)%vf, & - y_accel) - if (p > 0) then - call s_derive_acceleration_component(3, q_prim_ts(0)%vf, & - q_prim_ts(1)%vf, & - q_prim_ts(2)%vf, & - q_prim_ts(3)%vf, & - z_accel) - end if + z_accel) end if + end if - do k = 0, p - do j = 0, n - do i = 0, m - if (p > 0) then - accel_mag(i, j, k) = sqrt(x_accel(i, j, k)**2d0 + & - y_accel(i, j, k)**2d0 + & - z_accel(i, j, k)**2d0) - elseif (n > 0) then - accel_mag(i, j, k) = sqrt(x_accel(i, j, k)**2d0 + & - y_accel(i, j, k)**2d0) - else - accel_mag(i, j, k) = x_accel(i, j, k) - end if - end do + do k = 0, p + do j = 0, n + do i = 0, m + if (p > 0) then + accel_mag(i, j, k) = sqrt(x_accel(i, j, k)**2d0 + & + y_accel(i, j, k)**2d0 + & + z_accel(i, j, k)**2d0) + elseif (n > 0) then + accel_mag(i, j, k) = sqrt(x_accel(i, j, k)**2d0 + & + y_accel(i, j, k)**2d0) + else + accel_mag(i, j, k) = x_accel(i, j, k) + end if end do end do + end do - call s_write_probe_files(t_step, q_cons_ts(1)%vf, accel_mag) - end if + call s_write_probe_files(t_step, q_cons_ts(1)%vf, accel_mag) end if end subroutine s_compute_derived_variables ! --------------------------- @@ -373,1035 +348,12 @@ subroutine s_derive_acceleration_component(i, q_prim_vf, q_prim_vf1, & end subroutine s_derive_acceleration_component ! -------------------------- - !> This subroutine is used together with the volume fraction - !! model and when called upon, it computes the location of - !! of the center of mass for each fluid from the inputted - !! primitive variables, q_prim_vf. The computed location - !! is then written to a formatted data file by the root - !! process. - !! @param q_prim_vf Primitive variables - !! @param q_prim_vf1 Primitive variables - !! @param q_prim_vf2 Primitive variables - !! @param q_prim_vf3 Primitive variables - !! @param q_com Mass,x-location,y-location,z-location,x-velocity,y-velocity,z-velocity, - !! x-acceleration, y-acceleration, z-acceleration, weighted - subroutine s_derive_center_of_mass(q_prim_vf, q_prim_vf1, q_prim_vf2, q_prim_vf3, q_com) - - type(scalar_field), dimension(sys_size), intent(IN) :: q_prim_vf - type(scalar_field), dimension(sys_size), intent(IN) :: q_prim_vf1 - type(scalar_field), dimension(sys_size), intent(IN) :: q_prim_vf2 - type(scalar_field), dimension(sys_size), intent(IN) :: q_prim_vf3 - real(kind(0d0)), dimension(num_fluids, 10), intent(INOUT) :: q_com - - real(kind(0d0)) :: xbeg, xend, ybeg, yend, zbeg, zend !< - !! Maximum and minimum values of cell boundaries in each direction used in check for - !! reflective BC in computation of center of mass - - integer :: i, j, k, l !< Generic loop iterators - - real(kind(0d0)) :: tmp !< Temporary variable to store quantity for mpi_allreduce - - real(kind(0d0)) :: dV !< Discrete cell volume - - real(kind(0d0)) :: cart_u_x, cart_u_y, & - cart_u_x1, cart_u_y1, & - cart_u_x2, cart_u_y2, & - cart_u_x3, cart_u_y3 !< - !! Cartesian velocities - - if (n == 0) then !1D simulation - - do i = 1, num_fluids !Loop over individual fluids - if (com_wrt(i)) then - q_com(i, :) = 0d0 - do l = 0, p !Loop over grid - do k = 0, n - do j = 0, m - - dV = dx(j) - - ! Mass - q_com(i, 1) = q_com(i, 1) + q_prim_vf(i)%sf(j, k, l)*dV - ! x-location weighted - q_com(i, 2) = q_com(i, 2) + q_prim_vf(i)%sf(j, k, l)*dV*x_cc(j) - ! x-velocity weighted - q_com(i, 5) = q_com(i, 5) + q_prim_vf(i)%sf(j, k, l)*dV*q_prim_vf(mom_idx%beg)%sf(j, k, l) - ! x-acceleration weighted - q_com(i, 8) = q_com(i, 8) + dV*(11d0*(q_prim_vf(i)%sf(j, k, l) & - *q_prim_vf(mom_idx%beg)%sf(j, k, l)) & - - 18d0*(q_prim_vf1(i)%sf(j, k, l)*q_prim_vf1(mom_idx%beg)%sf(j, k, l)) & - + 9d0*(q_prim_vf2(i)%sf(j, k, l)*q_prim_vf2(mom_idx%beg)%sf(j, k, l)) & - - 2d0*(q_prim_vf3(i)%sf(j, k, l)*q_prim_vf3(mom_idx%beg)%sf(j, k, l)))/(6d0*dt) - end do - end do - end do - ! Sum all components across all processors using MPI_ALLREDUCE - if (num_procs > 1) then - tmp = q_com(i, 1) - call s_mpi_allreduce_sum(tmp, q_com(i, 1)) - tmp = q_com(i, 2) - call s_mpi_allreduce_sum(tmp, q_com(i, 2)) - tmp = q_com(i, 5) - call s_mpi_allreduce_sum(tmp, q_com(i, 5)) - tmp = q_com(i, 8) - call s_mpi_allreduce_sum(tmp, q_com(i, 8)) - end if - - ! Compute quotients - q_com(i, 2) = q_com(i, 2)/q_com(i, 1) - q_com(i, 5) = q_com(i, 5)/q_com(i, 1) - q_com(i, 8) = q_com(i, 8)/q_com(i, 1) - end if - end do - - elseif (p == 0) then !2D simulation - - do i = 1, num_fluids !Loop over individual fluids - if (com_wrt(i)) then - q_com(i, :) = 0d0 - do l = 0, p !Loop over grid - do k = 0, n - do j = 0, m - - dV = dx(j)*dy(k) - - ! Mass - q_com(i, 1) = q_com(i, 1) + q_prim_vf(i)%sf(j, k, l)*dV - ! x-location weighted - q_com(i, 2) = q_com(i, 2) + q_prim_vf(i)%sf(j, k, l)*dV*x_cc(j) - ! y-location weighted - q_com(i, 3) = q_com(i, 3) + q_prim_vf(i)%sf(j, k, l)*dV*y_cc(k) - ! x-velocity weighted - q_com(i, 5) = q_com(i, 5) + q_prim_vf(i)%sf(j, k, l)*dV*q_prim_vf(mom_idx%beg)%sf(j, k, l) - ! y-velocity weighted - q_com(i, 6) = q_com(i, 6) + q_prim_vf(i)%sf(j, k, l)*dV*q_prim_vf(mom_idx%beg + 1)%sf(j, k, l) - ! x-acceleration weighted - q_com(i, 8) = q_com(i, 8) + dV* & - (11d0*(q_prim_vf(i)%sf(j, k, l)*q_prim_vf(mom_idx%beg)%sf(j, k, l)) & - - 18d0*(q_prim_vf1(i)%sf(j, k, l)*q_prim_vf1(mom_idx%beg)%sf(j, k, l)) & - + 9d0*(q_prim_vf2(i)%sf(j, k, l)*q_prim_vf2(mom_idx%beg)%sf(j, k, l)) & - - 2d0*(q_prim_vf3(i)%sf(j, k, l)*q_prim_vf3(mom_idx%beg)%sf(j, k, l)))/(6d0*dt) - ! y-acceleration weighted - q_com(i, 9) = q_com(i, 9) + dV* & - (11d0*(q_prim_vf(i)%sf(j, k, l)*q_prim_vf(mom_idx%beg + 1)%sf(j, k, l)) & - - 18d0*(q_prim_vf1(i)%sf(j, k, l)*q_prim_vf1(mom_idx%beg + 1)%sf(j, k, l)) & - + 9d0*(q_prim_vf2(i)%sf(j, k, l)*q_prim_vf2(mom_idx%beg + 1)%sf(j, k, l)) & - - 2d0*(q_prim_vf3(i)%sf(j, k, l)*q_prim_vf3(mom_idx%beg + 1)%sf(j, k, l)))/(6d0*dt) - end do - end do - end do - ! Sum all components across all processors using MPI_ALLREDUCE - if (num_procs > 1) then - tmp = q_com(i, 1) - call s_mpi_allreduce_sum(tmp, q_com(i, 1)) - tmp = q_com(i, 2) - call s_mpi_allreduce_sum(tmp, q_com(i, 2)) - tmp = q_com(i, 3) - call s_mpi_allreduce_sum(tmp, q_com(i, 3)) - tmp = q_com(i, 5) - call s_mpi_allreduce_sum(tmp, q_com(i, 5)) - tmp = q_com(i, 6) - call s_mpi_allreduce_sum(tmp, q_com(i, 6)) - tmp = q_com(i, 8) - call s_mpi_allreduce_sum(tmp, q_com(i, 8)) - tmp = q_com(i, 9) - call s_mpi_allreduce_sum(tmp, q_com(i, 9)) - end if - - ! Compute quotients - q_com(i, 2) = q_com(i, 2)/q_com(i, 1) - q_com(i, 3) = q_com(i, 3)/q_com(i, 1) - q_com(i, 5) = q_com(i, 5)/q_com(i, 1) - q_com(i, 6) = q_com(i, 6)/q_com(i, 1) - q_com(i, 8) = q_com(i, 8)/q_com(i, 1) - q_com(i, 9) = q_com(i, 9)/q_com(i, 1) - end if - end do - - else !3D simulation - - do i = 1, num_fluids !Loop over individual fluids - if (com_wrt(i)) then - q_com(i, :) = 0d0 - do l = 0, p !Loop over grid - do k = 0, n - do j = 0, m - if (grid_geometry == 3) then - - dV = (2d0*y_cb(k - 1)*dy(k) + dy(k)**2d0)/2d0*dx(j)*dz(l) - cart_u_x = q_prim_vf(mom_idx%beg + 1)%sf(j, k, l)*cos(z_cc(l)) - & - q_prim_vf(mom_idx%end)%sf(j, k, l)*sin(z_cc(l)) - cart_u_y = q_prim_vf(mom_idx%beg + 1)%sf(j, k, l)*sin(z_cc(l)) + & - q_prim_vf(mom_idx%end)%sf(j, k, l)*cos(z_cc(l)) - cart_u_x1 = q_prim_vf1(mom_idx%beg + 1)%sf(j, k, l)*cos(z_cc(l)) - & - q_prim_vf1(mom_idx%end)%sf(j, k, l)*sin(z_cc(l)) - cart_u_y1 = q_prim_vf1(mom_idx%beg + 1)%sf(j, k, l)*sin(z_cc(l)) + & - q_prim_vf1(mom_idx%end)%sf(j, k, l)*cos(z_cc(l)) - cart_u_x2 = q_prim_vf2(mom_idx%beg + 1)%sf(j, k, l)*cos(z_cc(l)) - & - q_prim_vf2(mom_idx%end)%sf(j, k, l)*sin(z_cc(l)) - cart_u_y2 = q_prim_vf2(mom_idx%beg + 1)%sf(j, k, l)*sin(z_cc(l)) + & - q_prim_vf2(mom_idx%end)%sf(j, k, l)*cos(z_cc(l)) - cart_u_x3 = q_prim_vf3(mom_idx%beg + 1)%sf(j, k, l)*cos(z_cc(l)) - & - q_prim_vf3(mom_idx%end)%sf(j, k, l)*sin(z_cc(l)) - cart_u_y3 = q_prim_vf3(mom_idx%beg + 1)%sf(j, k, l)*sin(z_cc(l)) + & - q_prim_vf3(mom_idx%end)%sf(j, k, l)*cos(z_cc(l)) - - ! Mass - q_com(i, 1) = q_com(i, 1) + q_prim_vf(i)%sf(j, k, l)*dV - ! x-location weighted - q_com(i, 2) = q_com(i, 2) + q_prim_vf(i)%sf(j, k, l)*dV*y_cc(k)*cos(z_cc(l)) - ! y-location weighted - q_com(i, 3) = q_com(i, 3) + q_prim_vf(i)%sf(j, k, l)*dV*y_cc(k)*sin(z_cc(l)) - ! z-location weighted - q_com(i, 4) = q_com(i, 4) + q_prim_vf(i)%sf(j, k, l)*dV*x_cc(j) - ! x-velocity weighted - q_com(i, 5) = q_com(i, 5) + q_prim_vf(i)%sf(j, k, l)*dV*cart_u_x - ! y-velocity weighted - q_com(i, 6) = q_com(i, 6) + q_prim_vf(i)%sf(j, k, l)*dV*cart_u_y - ! z-velocity weighted - q_com(i, 7) = q_com(i, 7) + q_prim_vf(i)%sf(j, k, l)*dV*q_prim_vf(mom_idx%beg)%sf(j, k, l) - ! x-acceleration weighted - q_com(i, 8) = q_com(i, 8) + dV* & - (11d0*(q_prim_vf(i)%sf(j, k, l)*cart_u_x) & - - 18d0*(q_prim_vf1(i)%sf(j, k, l)*cart_u_x1) & - + 9d0*(q_prim_vf2(i)%sf(j, k, l)*cart_u_x2) & - - 2d0*(q_prim_vf3(i)%sf(j, k, l)*cart_u_x3))/(6d0*dt) - ! y-acceleration weighted - q_com(i, 9) = q_com(i, 9) + dV* & - (11d0*(q_prim_vf(i)%sf(j, k, l)*cart_u_y) & - - 18d0*(q_prim_vf1(i)%sf(j, k, l)*cart_u_y1) & - + 9d0*(q_prim_vf2(i)%sf(j, k, l)*cart_u_y2) & - - 2d0*(q_prim_vf3(i)%sf(j, k, l)*cart_u_y3))/(6d0*dt) - ! z-acceleration weighted - q_com(i, 10) = q_com(i, 10) + dV* & - (11d0*(q_prim_vf(i)%sf(j, k, l)*q_prim_vf(mom_idx%beg)%sf(j, k, l)) & - - 18d0*(q_prim_vf1(i)%sf(j, k, l)*q_prim_vf1(mom_idx%beg)%sf(j, k, l)) & - + 9d0*(q_prim_vf2(i)%sf(j, k, l)*q_prim_vf2(mom_idx%beg)%sf(j, k, l)) & - - 2d0*(q_prim_vf3(i)%sf(j, k, l)*q_prim_vf3(mom_idx%beg)%sf(j, k, l)))/(6d0*dt) - else - - dV = dx(j)*dy(k)*dz(l) - - ! Mass - q_com(i, 1) = q_com(i, 1) + q_prim_vf(i)%sf(j, k, l)*dV - ! x-location weighted - q_com(i, 2) = q_com(i, 2) + q_prim_vf(i)%sf(j, k, l)*dV*x_cc(j) - ! y-location weighted - q_com(i, 3) = q_com(i, 3) + q_prim_vf(i)%sf(j, k, l)*dV*y_cc(k) - ! z-location weighted - q_com(i, 4) = q_com(i, 4) + q_prim_vf(i)%sf(j, k, l)*dV*z_cc(l) - ! x-velocity weighted - q_com(i, 5) = q_com(i, 5) + q_prim_vf(i)%sf(j, k, l)*dV*q_prim_vf(mom_idx%beg)%sf(j, k, l) - ! y-velocity weighted - q_com(i, 6) = q_com(i, 6) + q_prim_vf(i)%sf(j, k, l)*dV*q_prim_vf(mom_idx%beg + 1)%sf(j, k, l) - ! z-velocity weighted - q_com(i, 7) = q_com(i, 7) + q_prim_vf(i)%sf(j, k, l)*dV*q_prim_vf(mom_idx%end)%sf(j, k, l) - ! x-acceleration weighted - q_com(i, 8) = q_com(i, 8) + dV* & - (11d0*(q_prim_vf(i)%sf(j, k, l)*q_prim_vf(mom_idx%beg)%sf(j, k, l)) & - - 18d0*(q_prim_vf1(i)%sf(j, k, l)*q_prim_vf1(mom_idx%beg)%sf(j, k, l)) & - + 9d0*(q_prim_vf2(i)%sf(j, k, l)*q_prim_vf2(mom_idx%beg)%sf(j, k, l)) & - - 2d0*(q_prim_vf3(i)%sf(j, k, l)*q_prim_vf3(mom_idx%beg)%sf(j, k, l)))/(6d0*dt) - ! y-acceleration weighted - q_com(i, 9) = q_com(i, 9) + dV* & - (11d0*(q_prim_vf(i)%sf(j, k, l)*q_prim_vf(mom_idx%beg + 1)%sf(j, k, l)) & - - 18d0*(q_prim_vf1(i)%sf(j, k, l)*q_prim_vf1(mom_idx%beg + 1)%sf(j, k, l)) & - + 9d0*(q_prim_vf2(i)%sf(j, k, l)*q_prim_vf2(mom_idx%beg + 1)%sf(j, k, l)) & - - 2d0*(q_prim_vf3(i)%sf(j, k, l)*q_prim_vf3(mom_idx%beg + 1)%sf(j, k, l)))/(6d0*dt) - ! z-acceleration weighted - q_com(i, 10) = q_com(i, 10) + dV* & - (11d0*(q_prim_vf(i)%sf(j, k, l)*q_prim_vf(mom_idx%end)%sf(j, k, l)) & - - 18d0*(q_prim_vf1(i)%sf(j, k, l)*q_prim_vf1(mom_idx%end)%sf(j, k, l)) & - + 9d0*(q_prim_vf2(i)%sf(j, k, l)*q_prim_vf2(mom_idx%end)%sf(j, k, l)) & - - 2d0*(q_prim_vf3(i)%sf(j, k, l)*q_prim_vf3(mom_idx%end)%sf(j, k, l)))/(6d0*dt) - end if - end do - end do - end do - ! Sum all components across all processors using MPI_ALLREDUCE - if (num_procs > 1) then - tmp = q_com(i, 1) - call s_mpi_allreduce_sum(tmp, q_com(i, 1)) - tmp = q_com(i, 2) - call s_mpi_allreduce_sum(tmp, q_com(i, 2)) - tmp = q_com(i, 3) - call s_mpi_allreduce_sum(tmp, q_com(i, 3)) - tmp = q_com(i, 4) - call s_mpi_allreduce_sum(tmp, q_com(i, 4)) - tmp = q_com(i, 5) - call s_mpi_allreduce_sum(tmp, q_com(i, 5)) - tmp = q_com(i, 6) - call s_mpi_allreduce_sum(tmp, q_com(i, 6)) - tmp = q_com(i, 7) - call s_mpi_allreduce_sum(tmp, q_com(i, 7)) - tmp = q_com(i, 8) - call s_mpi_allreduce_sum(tmp, q_com(i, 8)) - tmp = q_com(i, 9) - call s_mpi_allreduce_sum(tmp, q_com(i, 9)) - tmp = q_com(i, 10) - call s_mpi_allreduce_sum(tmp, q_com(i, 10)) - end if - - ! Compute quotients - q_com(i, 2) = q_com(i, 2)/q_com(i, 1) - q_com(i, 3) = q_com(i, 3)/q_com(i, 1) - q_com(i, 4) = q_com(i, 4)/q_com(i, 1) - q_com(i, 5) = q_com(i, 5)/q_com(i, 1) - q_com(i, 6) = q_com(i, 6)/q_com(i, 1) - q_com(i, 7) = q_com(i, 7)/q_com(i, 1) - q_com(i, 8) = q_com(i, 8)/q_com(i, 1) - q_com(i, 9) = q_com(i, 9)/q_com(i, 1) - q_com(i, 10) = q_com(i, 10)/q_com(i, 1) - end if - end do - - end if - - ! Find computational domain boundaries - if (num_procs > 1) then - call s_mpi_allreduce_min(minval(x_cb(-1:m)), xbeg) - call s_mpi_allreduce_max(maxval(x_cb(-1:m)), xend) - if (n > 0) then - call s_mpi_allreduce_min(minval(y_cb(-1:n)), ybeg) - call s_mpi_allreduce_max(maxval(y_cb(-1:n)), yend) - if (p > 0) then - call s_mpi_allreduce_min(minval(z_cb(-1:p)), zbeg) - call s_mpi_allreduce_max(maxval(z_cb(-1:p)), zend) - end if - end if - else - xbeg = minval(x_cb(-1:m)) - xend = maxval(x_cb(-1:m)) - if (n > 0) then - ybeg = minval(y_cb(-1:n)) - yend = maxval(y_cb(-1:n)) - if (p > 0) then - zbeg = minval(z_cb(-1:p)) - zend = maxval(z_cb(-1:p)) - end if - end if - end if - - do i = 1, num_fluids - if (com_wrt(i)) then - ! Check for reflective BC in x-direction - if (bc_x_glb%beg == -2) then - q_com(i, 1) = q_com(i, 1)*2d0 - q_com(i, 2) = xbeg - q_com(i, 5) = 0d0 - q_com(i, 8) = 0d0 - elseif (bc_x_glb%end == -2) then - q_com(i, 1) = q_com(i, 1)*2d0 - q_com(i, 2) = xend - q_com(i, 5) = 0d0 - q_com(i, 8) = 0d0 - end if - if (n > 0) then - ! Check for reflective BC in y-direction - if (bc_y_glb%beg == -2) then - q_com(i, 1) = q_com(i, 1)*2d0 - q_com(i, 3) = ybeg - q_com(i, 6) = 0d0 - q_com(i, 9) = 0d0 - elseif (bc_y_glb%end == -2) then - q_com(i, 1) = q_com(i, 1)*2d0 - q_com(i, 3) = yend - q_com(i, 6) = 0d0 - q_com(i, 9) = 0d0 - end if - if (p > 0) then - ! Check for reflective BC in z-direction - if (bc_z_glb%beg == -2) then - q_com(i, 1) = q_com(i, 1)*2d0 - q_com(i, 4) = zbeg - q_com(i, 7) = 0d0 - q_com(i, 10) = 0d0 - elseif (bc_z_glb%end == -2) then - q_com(i, 1) = q_com(i, 1)*2d0 - q_com(i, 4) = zend - q_com(i, 7) = 0d0 - q_com(i, 10) = 0d0 - end if - - end if - end if - end if - end do - - end subroutine s_derive_center_of_mass ! ---------------------------------- - - !> Subroutine to compute the higher moments in an attempt to find - !! the maximal size of the droplet - !! @param q_prim_vf Primitive variables - !! @param moments Higher moments (2 lateral directions, 5 moment orders) - subroutine s_derive_higher_moments(q_prim_vf, moments) - - type(scalar_field), dimension(sys_size), intent(IN) :: q_prim_vf - real(kind(0d0)), dimension(num_fluids, 2, 5), intent(INOUT) :: moments - - integer :: i, r !< Generic loop iterators - - ! Using the global boundary conditions, determine method of computing - ! higher moments for y-direction - if (n > 0) then - if ((bc_y_glb%beg /= -2) .and. (bc_y_glb%end /= -2)) then - ! Non-symmetric moments - call s_non_symmetric_moments(q_prim_vf, moments, 1) - elseif (((bc_y_glb%beg == -2) .and. (bc_y_glb%end == -2)) & - .or. & - ((bc_y_glb%beg == -1) .and. (bc_y_glb%end == -1))) then - print '(A)', 'Periodic boundary conditions in y-direction. '// & - 'Cannot compute higher moments. Exiting...' - call s_mpi_abort() - else - call s_symmetric_moments(q_prim_vf, moments, 1) - end if - - if (p > 0) then - if ((bc_z_glb%beg /= -2) .and. (bc_z_glb%end /= -2)) then - ! Non-symmetric moments - call s_non_symmetric_moments(q_prim_vf, moments, 2) - elseif (((bc_z_glb%beg == -2) .and. (bc_z_glb%end == -2)) & - .or. & - ((bc_z_glb%beg == -1) .and. (bc_z_glb%end == -1))) then - print '(A)', 'Periodic boundary conditions in z-direction. '// & - 'Cannot compute higher moments. Exiting...' - call s_mpi_abort() - else - call s_symmetric_moments(q_prim_vf, moments, 2) - end if - end if - else !1D simulation - do i = 1, num_fluids !Loop over individual fluids - if (com_wrt(i)) then - do r = 1, 5 - if (moment_order(r) /= dflt_int) then - moments(i, :, r) = 0d0 - else - moments(i, :, r) = dflt_real - end if - end do - end if - end do - end if - - end subroutine s_derive_higher_moments ! ----------------------------------------- - - !> Compute non-symmetric moments - !! @param q_prim_vf Primitive variables - !! @param moments Higher moments(2 lateral directions, 5 moment orders) - !! @param dir Current lateral direction - subroutine s_non_symmetric_moments(q_prim_vf, moments, dir) ! --------------------- - - type(scalar_field), dimension(sys_size), intent(IN) :: q_prim_vf - real(kind(0d0)), dimension(num_fluids, 2, 5), intent(INOUT) :: moments - integer, intent(IN) :: dir - - real(kind(0d0)), dimension(num_fluids, 5) :: pos_numer, neg_numer, pos_denom, neg_denom !< - !! Numerator and denominator place holders for computation - - real(kind(0d0)) :: numer_weight !< Numerator weight - real(kind(0d0)) :: main_term !< Constant term in both numerator and denominator - real(kind(0d0)) :: dV !< Discrete cell volume - real(kind(0d0)) :: cart_x, cart_y !< Cartesian x- and y-locations - integer :: i, j, k, l, r !< Generic loop iterators - real(kind(0d0)) :: tmp !< Temporary variable to store quantity for mpi_allreduce - - do i = 1, num_fluids - if (com_wrt(i)) then - pos_numer(i, :) = 0d0 - neg_numer(i, :) = 0d0 - pos_denom(i, :) = 0d0 - neg_denom(i, :) = 0d0 - do r = 1, 5 - if (moment_order(r) /= dflt_int) then - do l = 0, p - do k = 0, n - do j = 0, m - if (q_prim_vf(i)%sf(j, k, l) > 5d-1) then - if (grid_geometry == 3) then - dV = (2d0*y_cb(k - 1)*dy(k) + dy(k)**2d0)/2d0*dx(j)*dz(l) - cart_x = y_cc(k)*cos(z_cc(l)) - cart_y = y_cc(k)*sin(z_cc(l)) - main_term = q_prim_vf(i + E_idx)%sf(j, k, l)* & - (1d0 - q_prim_vf(i + E_idx)%sf(j, k, l))*dV - if ((dir == 1) .and. (cart_x >= 0d0)) then - numer_weight = cart_x**moment_order(r) - - pos_numer(i, r) = pos_numer(i, r) + numer_weight*main_term - pos_denom(i, r) = pos_denom(i, r) + main_term - elseif ((dir == 1) .and. (cart_x < 0d0)) then - numer_weight = cart_x**moment_order(r) - - neg_numer(i, r) = neg_numer(i, r) + numer_weight*main_term - neg_denom(i, r) = neg_denom(i, r) + main_term - elseif ((dir == 2) .and. (cart_y >= 0d0)) then - numer_weight = cart_y**moment_order(r) - - pos_numer(i, r) = pos_numer(i, r) + numer_weight*main_term - pos_denom(i, r) = pos_denom(i, r) + main_term - elseif ((dir == 2) .and. (cart_y < 0d0)) then - numer_weight = cart_y**moment_order(r) - - neg_numer(i, r) = neg_numer(i, r) + numer_weight*main_term - neg_denom(i, r) = neg_denom(i, r) + main_term - end if - else - if (n > 0) then - main_term = q_prim_vf(i + E_idx)%sf(j, k, l)* & - (1d0 - q_prim_vf(i + E_idx)%sf(j, k, l))* & - dx(j)*dy(k) - if (p > 0) then - main_term = main_term*dz(l) - end if - end if - if ((dir == 1) .and. (y_cc(k) >= 0d0)) then - numer_weight = y_cc(k)**moment_order(r) - - pos_numer(i, r) = pos_numer(i, r) + numer_weight*main_term - pos_denom(i, r) = pos_denom(i, r) + main_term - elseif ((dir == 1) .and. (y_cc(k) < 0d0)) then - numer_weight = y_cc(k)**moment_order(r) - - neg_numer(i, r) = neg_numer(i, r) + numer_weight*main_term - neg_denom(i, r) = neg_denom(i, r) + main_term - elseif ((dir == 2) .and. (z_cc(l) >= 0d0)) then - numer_weight = z_cc(l)**moment_order(r) - - pos_numer(i, r) = pos_numer(i, r) + numer_weight*main_term - pos_denom(i, r) = pos_denom(i, r) + main_term - elseif ((dir == 2) .and. (z_cc(l) < 0d0)) then - numer_weight = z_cc(l)**moment_order(r) - - neg_numer(i, r) = neg_numer(i, r) + numer_weight*main_term - neg_denom(i, r) = neg_denom(i, r) + main_term - end if - end if - end if - end do - end do - end do - ! Sum all components across all procs using MPI_ALLREDUCE - if (num_procs > 1) then - tmp = pos_numer(i, r) - call s_mpi_allreduce_sum(tmp, pos_numer(i, r)) - tmp = neg_numer(i, r) - call s_mpi_allreduce_sum(tmp, neg_numer(i, r)) - tmp = pos_denom(i, r) - call s_mpi_allreduce_sum(tmp, pos_denom(i, r)) - tmp = neg_denom(i, r) - call s_mpi_allreduce_sum(tmp, neg_denom(i, r)) - end if - ! Compute quotients and sum to get total moment - moments(i, dir, r) = (pos_numer(i, r)/pos_denom(i, r))**(1d0/moment_order(r)) + & - (neg_numer(i, r)/neg_denom(i, r))**(1d0/moment_order(r)) - else - moments(i, dir, r) = dflt_real - end if - end do - end if - end do - - end subroutine s_non_symmetric_moments ! ------------------------------------------ - - !> Compute symmetric moments - !! @param q_prim_vf Primitive variables - !! @param moments Higher moments(2 lateral directions, 5 moment orders) - !! @param dir Current lateral direction - subroutine s_symmetric_moments(q_prim_vf, moments, dir) ! --------------------- - - type(scalar_field), dimension(sys_size), intent(IN) :: q_prim_vf - real(kind(0d0)), dimension(num_fluids, 2, 5), intent(INOUT) :: moments - integer, intent(IN) :: dir - - real(kind(0d0)), dimension(num_fluids, 5) :: numer, denom !< - !! Numerator and denominator place holders for computation - - real(kind(0d0)) :: numer_weight !< Numerator weight - real(kind(0d0)) :: main_term !< Constant term in both numerator and denominator - real(kind(0d0)) :: dV !< Discrete cell volume - real(kind(0d0)) :: cart_x, cart_y !< Cartesian x- and y-locations - real(kind(0d0)) :: tmp !< Temporary variable to store quantity for mpi_allreduce - - integer :: i, j, k, l, r !< Generic loop iterators - - do i = 1, num_fluids - if (com_wrt(i)) then - numer(i, :) = 0d0 - denom(i, :) = 0d0 - do r = 1, 5 - if (moment_order(r) /= dflt_int) then - do l = 0, p - do k = 0, n - do j = 0, m - if (q_prim_vf(i)%sf(j, k, l) > 5d-1) then - if (grid_geometry == 3) then - dV = (2d0*y_cb(k - 1)*dy(k) + dy(k)**2d0)/2d0*dx(j)*dz(l) - cart_x = y_cc(k)*cos(z_cc(l)) - cart_y = y_cc(k)*sin(z_cc(l)) - main_term = q_prim_vf(i + E_idx)%sf(j, k, l)* & - (1d0 - q_prim_vf(i + E_idx)%sf(j, k, l))*dV - if (dir == 1) then - numer_weight = cart_x**moment_order(r) - - numer(i, r) = numer(i, r) + numer_weight*main_term - denom(i, r) = denom(i, r) + main_term - elseif (dir == 2) then - numer_weight = cart_y**moment_order(r) - - numer(i, r) = numer(i, r) + numer_weight*main_term - denom(i, r) = denom(i, r) + main_term - end if - else - if (n > 0) then - main_term = q_prim_vf(i + E_idx)%sf(j, k, l)* & - (1d0 - q_prim_vf(i + E_idx)%sf(j, k, l))* & - dx(j)*dy(k) - if (p > 0) then - main_term = main_term*dz(l) - end if - end if - if (dir == 1) then - numer_weight = y_cc(k)**moment_order(r) - - numer(i, r) = numer(i, r) + numer_weight*main_term - denom(i, r) = denom(i, r) + main_term - elseif (dir == 2) then - numer_weight = z_cc(l)**moment_order(r) - - numer(i, r) = numer(i, r) + numer_weight*main_term - denom(i, r) = denom(i, r) + main_term - end if - end if - end if - end do - end do - end do - ! Sum all components across all procs using MPI_ALLREDUCE - if (num_procs > 1) then - tmp = numer(i, r) - call s_mpi_allreduce_sum(tmp, numer(i, r)) - tmp = denom(i, r) - call s_mpi_allreduce_sum(tmp, denom(i, r)) - end if - ! Compute quotients and sum to get total moment - moments(i, dir, r) = (numer(i, r)/denom(i, r))**(1d0/moment_order(r)) - else - moments(i, dir, r) = dflt_real - end if - end do - end if - end do - - end subroutine s_symmetric_moments ! ------------------------------------------ - - !> This subroutine is used together with the volume fraction model - !! and when called upon, it computes the min and max bounds of the - !! fluid in each direction in the domain. - !! @param q_prim_vf Primitive variables - !! @param bounds Variables storing the min and max bounds of the fluids - subroutine s_derive_fluid_bounds(q_prim_vf, bounds) ! ----------------------- - - type(scalar_field), dimension(sys_size), intent(IN) :: q_prim_vf - real(kind(0d0)), dimension(num_fluids, 5, 6), intent(INOUT) :: bounds - - real(kind(0d0)) :: cart_x, cart_y, cart_z !< Cartesian x,y,z-locations - real(kind(0d0)) :: tmp !< Temporary variable to store quantity for mpi_allreduce - - integer :: i, j, k, l, r !< Generic loop iterators - - if (n == 0) then ! 1D simulation - do i = 1, num_fluids !Loop over individual fluids - if (cb_wrt(i)) then - bounds(i, :, 1) = -1d0*dflt_real ! 1d6 - bounds(i, :, 2) = dflt_real ! -1d6 - do r = 1, 5 - if (threshold_mf(r) /= dflt_real) then - do l = 0, p !Loop over grid - do k = 0, n - do j = 0, m - if ((q_prim_vf(i + E_idx)%sf(j, k, l) >= threshold_mf(r)) & - .and. (x_cb(j - 1) <= bounds(i, r, 1))) then - bounds(i, r, 1) = x_cb(j - 1) - elseif ((q_prim_vf(i + E_idx)%sf(j, k, l) >= threshold_mf(r)) & - .and. (x_cb(j) >= bounds(i, r, 2))) then - bounds(i, r, 2) = x_cb(j) - end if - end do - end do - end do - - if (num_procs > 1) then - tmp = bounds(i, r, 1) - call s_mpi_allreduce_min(tmp, bounds(i, r, 1)) - tmp = bounds(i, r, 2) - call s_mpi_allreduce_max(tmp, bounds(i, r, 2)) - end if - else - bounds(i, r, 1) = dflt_real - bounds(i, r, 2) = dflt_real - end if - end do - end if - end do - elseif (p == 0) then ! 2D simulation - do i = 1, num_fluids !Loop over individual fluids - if (cb_wrt(i)) then - bounds(i, :, 1) = -1d0*dflt_real ! 1d6 - bounds(i, :, 2) = dflt_real ! -1d6 - bounds(i, :, 3) = -1d0*dflt_real ! 1d6 - bounds(i, :, 4) = dflt_real ! -1d6 - do r = 1, 5 - if (threshold_mf(r) /= dflt_real) then - do l = 0, p ! Loop over grid - do k = 0, n - do j = 0, m - if ((q_prim_vf(i + E_idx)%sf(j, k, l) >= threshold_mf(r)) & - .and. (x_cb(j - 1) <= bounds(i, r, 1))) then - bounds(i, r, 1) = x_cb(j - 1) - elseif ((q_prim_vf(i + E_idx)%sf(j, k, l) >= threshold_mf(r)) & - .and. (x_cb(j) >= bounds(i, r, 2))) then - bounds(i, r, 2) = x_cb(j) - end if - if ((q_prim_vf(i + E_idx)%sf(j, k, l) >= threshold_mf(r)) & - .and. (y_cb(k - 1) <= bounds(i, r, 3))) then - bounds(i, r, 3) = y_cb(k - 1) - elseif ((q_prim_vf(i + E_idx)%sf(j, k, l) >= threshold_mf(r)) & - .and. (y_cb(k) >= bounds(i, r, 4))) then - bounds(i, r, 4) = y_cb(k) - end if - end do - end do - end do - - if (num_procs > 1) then - tmp = bounds(i, r, 1) - call s_mpi_allreduce_min(tmp, bounds(i, r, 1)) - tmp = bounds(i, r, 2) - call s_mpi_allreduce_max(tmp, bounds(i, r, 2)) - tmp = bounds(i, r, 3) - call s_mpi_allreduce_min(tmp, bounds(i, r, 3)) - tmp = bounds(i, r, 4) - call s_mpi_allreduce_max(tmp, bounds(i, r, 4)) - end if - else - bounds(i, r, 1) = dflt_real - bounds(i, r, 2) = dflt_real - bounds(i, r, 3) = dflt_real - bounds(i, r, 4) = dflt_real - end if - end do - end if - end do - else ! 3D simulation - do i = 1, num_fluids !Loop over individual fluids - if (cb_wrt(i)) then - bounds(i, :, 1) = -1d0*dflt_real ! 1d6 - bounds(i, :, 2) = dflt_real ! -1d6 - bounds(i, :, 3) = -1d0*dflt_real ! 1d6 - bounds(i, :, 4) = dflt_real ! -1d6 - bounds(i, :, 5) = -1d0*dflt_real ! 1d6 - bounds(i, :, 6) = dflt_real ! -1d6 - do r = 1, 5 - if (threshold_mf(r) /= dflt_real) then - do l = 0, p ! Loop over grid - do k = 0, n - do j = 0, m - if (grid_geometry == 3) then - cart_x = y_cc(k)*cos(z_cc(l)) - cart_y = y_cc(k)*sin(z_cc(l)) - cart_z = x_cc(j) - if ((q_prim_vf(i + E_idx)%sf(j, k, l) >= threshold_mf(r)) & - .and. (cart_x <= bounds(i, r, 1))) then - bounds(i, r, 1) = cart_x - elseif ((q_prim_vf(i + E_idx)%sf(j, k, l) >= threshold_mf(r)) & - .and. (cart_x >= bounds(i, r, 2))) then - bounds(i, r, 2) = cart_x - end if - if ((q_prim_vf(i + E_idx)%sf(j, k, l) >= threshold_mf(r)) & - .and. (cart_y <= bounds(i, r, 3))) then - bounds(i, r, 3) = cart_y - elseif ((q_prim_vf(i + E_idx)%sf(j, k, l) >= threshold_mf(r)) & - .and. (cart_y >= bounds(i, r, 4))) then - bounds(i, r, 4) = cart_y - end if - if ((q_prim_vf(i + E_idx)%sf(j, k, l) >= threshold_mf(r)) & - .and. (cart_z <= bounds(i, r, 5))) then - bounds(i, r, 5) = cart_z - elseif ((q_prim_vf(i + E_idx)%sf(j, k, l) >= threshold_mf(r)) & - .and. (cart_z >= bounds(i, r, 6))) then - bounds(i, r, 6) = cart_z - end if - else - if ((q_prim_vf(i + E_idx)%sf(j, k, l) >= threshold_mf(r)) & - .and. (x_cb(j - 1) <= bounds(i, r, 1))) then - bounds(i, r, 1) = x_cb(j - 1) - elseif ((q_prim_vf(i + E_idx)%sf(j, k, l) >= threshold_mf(r)) & - .and. (x_cb(j) >= bounds(i, r, 2))) then - bounds(i, r, 2) = x_cb(j) - end if - if ((q_prim_vf(i + E_idx)%sf(j, k, l) >= threshold_mf(r)) & - .and. (y_cb(k - 1) <= bounds(i, r, 3))) then - bounds(i, r, 3) = y_cb(k - 1) - elseif ((q_prim_vf(i + E_idx)%sf(j, k, l) >= threshold_mf(r)) & - .and. (y_cb(k) >= bounds(i, r, 4))) then - bounds(i, r, 4) = y_cb(k) - end if - if ((q_prim_vf(i + E_idx)%sf(j, k, l) >= threshold_mf(r)) & - .and. (z_cb(l - 1) <= bounds(i, r, 5))) then - bounds(i, r, 5) = z_cb(l - 1) - elseif ((q_prim_vf(i + E_idx)%sf(j, k, l) >= threshold_mf(r)) & - .and. (z_cb(l) >= bounds(i, r, 6))) then - bounds(i, r, 6) = z_cb(l) - end if - end if - end do - end do - end do - - if (num_procs > 1) then - tmp = bounds(i, r, 1) - call s_mpi_allreduce_min(tmp, bounds(i, r, 1)) - tmp = bounds(i, r, 2) - call s_mpi_allreduce_max(tmp, bounds(i, r, 2)) - tmp = bounds(i, r, 3) - call s_mpi_allreduce_min(tmp, bounds(i, r, 3)) - tmp = bounds(i, r, 4) - call s_mpi_allreduce_max(tmp, bounds(i, r, 4)) - tmp = bounds(i, r, 5) - call s_mpi_allreduce_min(tmp, bounds(i, r, 5)) - tmp = bounds(i, r, 6) - call s_mpi_allreduce_max(tmp, bounds(i, r, 6)) - end if - else - bounds(i, r, 1) = dflt_real - bounds(i, r, 2) = dflt_real - bounds(i, r, 3) = dflt_real - bounds(i, r, 4) = dflt_real - bounds(i, r, 5) = dflt_real - bounds(i, r, 6) = dflt_real - end if - end do - end if - end do - end if - - end subroutine s_derive_fluid_bounds ! ------------------------------------------- - - !> This subroutine is used together with the volume fraction model - !! and when called upon, it computes the total mass of a fluid in the - !! entire domain for which the volume fraction is greater than a - !! threshold value. This gives the mass of the coherent body. - !! @param q_prim_vf Primitive variables - !! @param cb_mass Coherent body mass - subroutine s_derive_coherent_body(q_prim_vf, cb_mass) ! -------------------------- - - type(scalar_field), dimension(sys_size), intent(IN) :: q_prim_vf - real(kind(0d0)), dimension(num_fluids, 10), intent(INOUT) :: cb_mass - - real(kind(0d0)) :: dV !< Discrete cell volume - integer :: i, j, k, l, r !< Generic loop iterators - real(kind(0d0)) :: tmp !< Temporary variable to store quantity for mpi_allreduce - - do i = 1, num_fluids !Loop over individual fluids - if (cb_wrt(i)) then - cb_mass(i, :) = 0d0 - do r = 1, 5 ! Volume fraction threshold values - if (threshold_mf(r) /= dflt_real) then - do l = 0, p !Loop over grid - do k = 0, n - do j = 0, m - if (q_prim_vf(i + E_idx)%sf(j, k, l) >= threshold_mf(r)) then - if (n == 0) then - dV = dx(j) - elseif (p == 0) then - dV = dx(j)*dy(k) - else - if (grid_geometry == 3) then - dV = (2d0*y_cb(k - 1)*dy(k) + dy(k)**2d0)/2d0*dx(j)*dz(l) - else - dV = dx(j)*dy(k)*dz(l) - end if - end if - cb_mass(i, r) = cb_mass(i, r) + q_prim_vf(i)%sf(j, k, l)*dV ! Mass - cb_mass(i, r + 5) = cb_mass(i, r + 5) + dV ! Volume - end if - end do - end do - end do - - if (num_procs > 1) then - tmp = cb_mass(i, r) - call s_mpi_allreduce_sum(tmp, cb_mass(i, r)) - tmp = cb_mass(i, r + 5) - call s_mpi_allreduce_sum(tmp, cb_mass(i, r + 5)) - end if - else - cb_mass(i, r) = dflt_real - cb_mass(i, r + 5) = dflt_real - end if - end do - end if - end do - - do i = 1, num_fluids - if (cb_wrt(i)) then - do r = 1, 5 - if (threshold_mf(r) /= dflt_real) then - ! Check for reflective BC in x-direction - if (bc_x_glb%beg == -2) then - cb_mass(i, r) = cb_mass(i, r)*2d0 - cb_mass(i, r + 5) = cb_mass(i, r + 5)*2d0 - elseif (bc_x_glb%end == -2) then - cb_mass(i, r) = cb_mass(i, r)*2d0 - cb_mass(i, r + 5) = cb_mass(i, r + 5)*2d0 - end if - if (n > 0) then - ! Check for reflective BC in y-direction - if (bc_y_glb%beg == -2) then - cb_mass(i, r) = cb_mass(i, r)*2d0 - cb_mass(i, r + 5) = cb_mass(i, r + 5)*2d0 - elseif (bc_y_glb%end == -2) then - cb_mass(i, r) = cb_mass(i, r)*2d0 - cb_mass(i, r + 5) = cb_mass(i, r + 5)*2d0 - end if - if (p > 0) then - ! Check for reflective BC in z-direction - if (bc_z_glb%beg == -2) then - cb_mass(i, r) = cb_mass(i, r)*2d0 - cb_mass(i, r + 5) = cb_mass(i, r + 5)*2d0 - elseif (bc_z_glb%end == -2) then - cb_mass(i, r) = cb_mass(i, r)*2d0 - cb_mass(i, r + 5) = cb_mass(i, r + 5)*2d0 - end if - - end if - end if - end if - end do - end if - end do - - end subroutine s_derive_coherent_body ! ------------------------------ - - !> This subroutine is used together with the volume fraction model - !! and when called upon, it computes the centerline length of the - !! fluid. - !! @param q_prim_vf Primitive variables - !! @param cntrline Variables storing the centerline length of the fluids - subroutine s_derive_centerline(q_prim_vf, cntrline) ! -------------------------- - - type(scalar_field), dimension(sys_size), intent(IN) :: q_prim_vf - real(kind(0d0)), dimension(num_fluids, 5), intent(INOUT) :: cntrline - - real(kind(0d0)), dimension(5) :: cntrmin, cntrmax !< Placeholders - real(kind(0d0)) :: tmp !< Temporary variable to store quantity for mpi_allreduce - - integer :: i, j, k, l, r !< Generic loop iterators - - if (n == 0) then ! 1D simulation - do i = 1, num_fluids - if (cb_wrt(i)) then - do r = 1, 5 - if (threshold_mf(r) /= dflt_real) then - cntrline(i, r) = 0d0 - else - cntrline(i, r) = dflt_real - end if - end do - end if - end do - elseif ((p == 0) .or. (grid_geometry == 3)) then ! 2D simulation - do i = 1, num_fluids - if (cb_wrt(i)) then - cntrmin(:) = -1d0*dflt_real - cntrmax(:) = dflt_real - do r = 1, 5 - if (threshold_mf(r) /= dflt_real) then - do l = 0, p - do k = 0, n - do j = 0, m - if ((y_cb(k - 1) == 0d0) .and. & - (q_prim_vf(i + E_idx)%sf(j, k, l) >= threshold_mf(r)) .and. & - (x_cb(j - 1) <= cntrmin(r))) then - cntrmin(r) = x_cb(j - 1) - elseif ((y_cb(k - 1) == 0d0) .and. & - (q_prim_vf(i + E_idx)%sf(j, k, l) >= threshold_mf(r)) .and. & - (x_cb(j) >= cntrmax(r))) then - cntrmax(r) = x_cb(j) - end if - end do - end do - end do - - if (num_procs > 1) then - tmp = cntrmin(r) - call s_mpi_allreduce_min(tmp, cntrmin(r)) - tmp = cntrmax(r) - call s_mpi_allreduce_max(tmp, cntrmax(r)) - end if - - cntrline(i, r) = cntrmax(r) - cntrmin(r) - else - cntrline(i, r) = dflt_real - end if - end do - end if - end do - else ! 3D simulation - do i = 1, num_fluids - if (cb_wrt(i)) then - cntrmin(:) = -1d0*dflt_real - cntrmax(:) = dflt_real - do r = 1, 5 - if (threshold_mf(r) /= dflt_real) then - do l = 0, p - do k = 0, n - do j = 0, m - if ((y_cb(k - 1) == 0d0) .and. & - (z_cb(l - 1) == 0d0) .and. & - (q_prim_vf(i + E_idx)%sf(j, k, l) >= threshold_mf(r)) .and. & - (x_cb(j - 1) <= cntrmin(r))) then - cntrmin(r) = x_cb(j - 1) - elseif ((y_cb(k - 1) == 0d0) .and. & - (z_cb(l - 1) == 0d0) .and. & - (q_prim_vf(i + E_idx)%sf(j, k, l) >= threshold_mf(r)) .and. & - (x_cb(j) >= cntrmax(r))) then - cntrmax(r) = x_cb(j) - end if - end do - end do - end do - - if (num_procs > 1) then - tmp = cntrmin(r) - call s_mpi_allreduce_min(tmp, cntrmin(r)) - tmp = cntrmax(r) - call s_mpi_allreduce_max(tmp, cntrmax(r)) - end if - - cntrline(i, r) = cntrmax(r) - cntrmin(r) - else - cntrline(i, r) = dflt_real - end if - end do - end if - end do - end if - - end subroutine s_derive_centerline ! ---------------------------------- !> Deallocation procedures for the module subroutine s_finalize_derived_variables_module() ! ------------------- ! Closing CoM and flow probe files if (proc_rank == 0) then - if (any(com_wrt)) then - call s_close_com_files() - end if - if (any(cb_wrt)) then - call s_close_cb_files() - end if if (probe_wrt) then call s_close_probe_files() end if diff --git a/src/simulation/m_global_parameters.fpp b/src/simulation/m_global_parameters.fpp index 054d45688c..71b42b1ed3 100644 --- a/src/simulation/m_global_parameters.fpp +++ b/src/simulation/m_global_parameters.fpp @@ -225,15 +225,12 @@ module m_global_parameters !! it is a measure of the half-size of the finite-difference stencil for the !! selected order of accuracy. - logical, dimension(num_fluids_max) :: com_wrt, cb_wrt logical :: probe_wrt logical :: integral_wrt integer :: num_probes integer :: num_integrals type(probe_parameters), dimension(num_probes_max) :: probe type(integral_parameters), dimension(num_probes_max) :: integral - real(kind(0d0)), dimension(5) :: threshold_mf - integer, dimension(5) :: moment_order !> @name Reference density and pressure for Tait EOS !> @{ @@ -444,8 +441,6 @@ contains end do fd_order = dflt_int - com_wrt = .false. - cb_wrt = .false. probe_wrt = .false. integral_wrt = .false. num_probes = dflt_int @@ -466,11 +461,6 @@ contains integral(i)%ymax = dflt_real end do - do i = 1, 5 - threshold_mf(i) = dflt_real - moment_order(i) = dflt_int - end do - end subroutine s_assign_default_values_to_user_inputs ! ---------------- !> The computation of parameters, the allocation of memory, diff --git a/src/simulation/m_mpi_proxy.fpp b/src/simulation/m_mpi_proxy.fpp index b5b85c5429..a543d74466 100644 --- a/src/simulation/m_mpi_proxy.fpp +++ b/src/simulation/m_mpi_proxy.fpp @@ -202,7 +202,6 @@ contains & 'bubble_model', 'thermal', 'R0_type', 'num_mono' ] call MPI_BCAST(${VAR}$, 1, MPI_INTEGER, 0, MPI_COMM_WORLD, ierr) #:endfor - call MPI_BCAST(moment_order(1), 5, MPI_INTEGER, 0, MPI_COMM_WORLD, ierr) #:for VAR in [ 'run_time_info','cyl_coord', 'adv_alphan', 'mpp_lim', & & 'mapped_weno', 'mp_weno', 'cu_mpi', 'weno_flat', 'riemann_flat', & @@ -211,14 +210,11 @@ contains & 'polydisperse', 'qbmm', 'monopole', 'probe_wrt', 'integral_wrt' ] call MPI_BCAST(${VAR}$, 1, MPI_LOGICAL, 0, MPI_COMM_WORLD, ierr) #:endfor - call MPI_BCAST(com_wrt(1), num_fluids_max, MPI_LOGICAL, 0, MPI_COMM_WORLD, ierr) - call MPI_BCAST(cb_wrt(1), num_fluids_max, MPI_LOGICAL, 0, MPI_COMM_WORLD, ierr) #:for VAR in [ 'dt','weno_eps','pref','rhoref','R0ref','Web','Ca', & & 'Re_inv','poly_sigma' ] call MPI_BCAST(${VAR}$, 1, MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr) #:endfor - call MPI_BCAST(threshold_mf(1), 5, MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr) #:if not MFC_CASE_OPTIMIZATION call MPI_BCAST(weno_order, 1, MPI_INTEGER, 0, MPI_COMM_WORLD, ierr) diff --git a/src/simulation/m_rhs.fpp b/src/simulation/m_rhs.fpp index 8c77940c78..62107e8530 100644 --- a/src/simulation/m_rhs.fpp +++ b/src/simulation/m_rhs.fpp @@ -1739,7 +1739,7 @@ contains end do ! END: Dimensional Splitting Loop ================================= - if (run_time_info .or. probe_wrt .or. any(com_wrt) .or. any(cb_wrt)) then + if (run_time_info .or. probe_wrt) then ix%beg = -buff_size; iy%beg = 0; iz%beg = 0 if (n > 0) iy%beg = -buff_size; diff --git a/src/simulation/m_start_up.fpp b/src/simulation/m_start_up.fpp index b03c0e628f..759edd69f0 100644 --- a/src/simulation/m_start_up.fpp +++ b/src/simulation/m_start_up.fpp @@ -82,9 +82,8 @@ contains riemann_solver, wave_speeds, avg_state, & bc_x, bc_y, bc_z, & hypoelasticity, & - fluid_pp, com_wrt, cb_wrt, probe_wrt, & + fluid_pp, probe_wrt, & fd_order, probe, num_probes, t_step_old, & - threshold_mf, moment_order, & alt_soundspeed, mixture_err, weno_Re_flux, & null_weights, precision, parallel_io, cyl_coord, & rhoref, pref, bubbles, bubble_model, & @@ -433,22 +432,6 @@ contains print '(A)', 'Unsupported combination of values of '// & 'bc_z%beg and bc_z%end. Exiting ...' call s_mpi_abort() - elseif ((any(threshold_mf /= dflt_real)) & - .and. & - (all(cb_wrt .neqv. .true.))) then - print '(A)', 'Unsupported combination of cb_wrt '// & - 'and threshold_mf. Exiting ...' - call s_mpi_abort() - elseif ((any(moment_order /= dflt_int)) & - .and. & - (all(com_wrt .neqv. .true.))) then - print '(A)', 'Unsupported combination of com_wrt '// & - 'and moment_order. Exiting ...' - call s_mpi_abort() - elseif (any(cb_wrt) .and. (all(threshold_mf == dflt_real))) then - print '(A)', 'Unsupported combination of cb_wrt '// & - 'and threshold_mf. Exiting ...' - call s_mpi_abort() elseif (model_eqns == 1 .and. alt_soundspeed) then print '(A)', 'Unsupported combination of model_eqns '// & 'and alt_soundspeed. Exiting ...' diff --git a/src/simulation/m_time_steppers.fpp b/src/simulation/m_time_steppers.fpp index 5399053166..873c51e66f 100644 --- a/src/simulation/m_time_steppers.fpp +++ b/src/simulation/m_time_steppers.fpp @@ -104,7 +104,7 @@ contains end do ! Allocating the cell-average primitive ts variables - if (any(com_wrt) .or. any(cb_wrt) .or. probe_wrt) then + if (probe_wrt) then @:ALLOCATE(q_prim_ts(0:3)) do i = 0, 3 @@ -198,7 +198,7 @@ contains print *, 'wrote runtime info' #endif - if (any(com_wrt) .or. any(cb_wrt) .or. probe_wrt) then + if (probe_wrt) then call s_time_step_cycling(t_step) end if @@ -262,7 +262,7 @@ contains call s_write_run_time_information(q_prim_vf, t_step) end if - if (any(com_wrt) .or. any(cb_wrt) .or. probe_wrt) then + if (probe_wrt) then call s_time_step_cycling(t_step) end if @@ -345,7 +345,7 @@ contains call s_write_run_time_information(q_prim_vf, t_step) end if - if (any(com_wrt) .or. any(cb_wrt) .or. probe_wrt) then + if (probe_wrt) then call s_time_step_cycling(t_step) end if @@ -490,7 +490,7 @@ contains @:DEALLOCATE(q_cons_ts) ! Deallocating the cell-average primitive ts variables - if (any(com_wrt) .or. any(cb_wrt) .or. probe_wrt) then + if (probe_wrt) then do i = 0, 3 do j = 1, sys_size @:DEALLOCATE(q_prim_ts(i)%vf(j)%sf) diff --git a/toolchain/mfc/run/case_dicts.py b/toolchain/mfc/run/case_dicts.py index f36872ac99..cbd0a6679f 100644 --- a/toolchain/mfc/run/case_dicts.py +++ b/toolchain/mfc/run/case_dicts.py @@ -70,10 +70,9 @@ 'riemann_solver', 'wave_speeds', 'avg_state', 'commute_err', 'split_err', 'alt_crv', 'alt_soundspeed', 'regularization', 'reg_eps', 'null_weights', 'mixture_err', 'tvd_riemann_flux', 'tvd_rhs_flux', 'tvd_wave_speeds', - 'flux_lim', 'We_riemann_flux', 'We_rhs_flux', 'We_src', 'We_wave_speeds', - 'lsq_deriv', 'fd_order', 'com_wrt', 'num_probes', 'probe_wrt', 'cb_wrt', - 'threshold_mf', 'moment_order', 'bubble_model', 'Monopole', 'num_mono', - 'qbmm', 'R0_type', 'integral_wrt', 'num_integrals', 'cu_mpi' + 'flux_lim', 'lsq_deriv', 'fd_order', 'num_probes', 'probe_wrt', + 'bubble_model', 'Monopole', 'num_mono', 'qbmm', 'R0_type', 'integral_wrt', + 'num_integrals', 'cu_mpi' ] for cmp in ["x", "y", "z"]: @@ -81,9 +80,6 @@ SIMULATION.append(f'bc_{cmp}%end') for wrt_id in range(1,10+1): - SIMULATION.append(f'com_wrt({wrt_id})') - SIMULATION.append(f'cb_wrt({wrt_id})') - for cmp in ["x", "y", "z"]: SIMULATION.append(f'probe_wrt({wrt_id})%{cmp}') @@ -91,12 +87,6 @@ for cmp in ["x", "y", "z"]: SIMULATION.append(f'probe({probe_id})%{cmp}') -for mf_id in range(1,5+1): - SIMULATION.append(f'threshold_mf({mf_id})') - -for order_id in range(1,5+1): - SIMULATION.append(f'moment_order({order_id})') - for f_id in range(1,10+1): for attribute in ["gamma", "pi_inf", "mul0", "ss", "pv", "gamma_v", "M_v", "mu_v", "k_v", "G"]: