diff --git a/docs/documentation/case.md b/docs/documentation/case.md index 5f1d739530..9e0536f913 100644 --- a/docs/documentation/case.md +++ b/docs/documentation/case.md @@ -482,9 +482,9 @@ It is recommended to set `weno_eps` to $10^{-6}$ for WENO-JS, and to $10^{-40}$ - `mp_weno` activates monotonicity preservation in the WENO reconstruction (MPWENO) such that the values of reconstructed variables do not reside outside the range spanned by WENO stencil ([Balsara and Shu, 2000](references.md); [Suresh and Huynh, 1997](references.md)). -- `muscl_order` specifies the order of the MUSCL scheme that is used for spatial reconstruction of variables by an integer of 1, or 2, that corresponds to the 1st, and 2nd order respectively. When using `muscl_order = 2`, `muscl_lim` must be defined. +- `muscl_order` specifies the order of the MUSCL scheme that is used for spatial reconstruction of variables by an integer of 1, or 2, that corresponds to the 1st, and 2nd order respectively. When using `muscl_order = 2`, `muscl_lim` must be defined. -- `muscl_lim` specifies the slope limiter that is used in 2nd order MUSCL Reconstruction by an integer from 1 through 5. +- `muscl_lim` specifies the slope limiter that is used in 2nd order MUSCL Reconstruction by an integer from 1 through 5. `muscl_lim = 1`, `2`, `3`, `4`, and `5` correspond to minmod, monotonized central, Van Albada, Van Leer, and SUPERBEE, respectively. - `int_comp` activates interface compression using THINC used in MUSCL Reconstruction, with control parameters (`ic_eps`, and `ic_beta`). @@ -599,6 +599,24 @@ To restart the simulation from $k$-th time step, see [Restarting Cases](running. | `output_partial_domain` | Logical | Output part of the domain | | `[x,y,z]_output%beg` | Real | Beginning of the output domain in the [x,y,z]-direction | | `[x,y,z]_output%end` | Real | End of the output domain in the [x,y,z]-direction | +| `lag_txt_wrt` | Logical | Write Lagrangian bubble data to `.dat` files | +| `lag_header` | Logical | Write header to Lagrangian bubble `.dat` files | +| `lag_db_wrt` | Logical | Write Lagrangian bubble data to silo/hdf5 database files | +| `lag_id_wrt` | Logical | Add the global bubble idea to the database file | +| `lag_pos_wrt` | Logical | Add the bubble position to the database file | +| `lag_pos_prev_wrt` | Logical | Add the previous bubble position to the database file | +| `lag_vel_wrt` | Logical | Add the bubble translational velocity to the database file | +| `lag_rad_wrt` | Logical | Add the bubble radius to the database file | +| `lag_rvel_wrt` | Logical | Add the bubble radial velocity to the database file | +| `lag_r0_wrt` | Logical | Add the bubble initial radius to the database file | +| `lag_rmax_wrt` | Logical | Add the bubble maximum radius to the database file | +| `lag_rmin_wrt` | Logical | Add the bubble minimum radius to the database file | +| `lag_dphidt_wrt` | Logical | Add the bubble subgrid velocity potential to the database file | +| `lag_pres_wrt` | Logical | Add the bubble pressure to the database file | +| `lag_mv_wrt` | Logical | Add the bubble vapor mass to the database file | +| `lag_mg_wrt` | Logical | Add the bubble gas mass to the database file | +| `lag_betaT_wrt` | Logical | Add the bubble heat flux model coefficient to the database file | +| `lag_betaC_wrt` | Logical | Add the bubble mass flux model coefficient to the database file | The table lists formatted database output parameters. The parameters define variables that are outputted from simulation and file types and formats of data as well as options for post-processing. @@ -628,7 +646,7 @@ If `file_per_process` is true, then pre_process, simulation, and post_process mu - `output_partial_domain` activates the output of part of the domain specified by `[x,y,z]_output%beg` and `[x,y,z]_output%end`. This is useful for large domains where only a portion of the domain is of interest. -It is not supported when `precision = 1` and `format = 1`. +It is not supported when `precision = 1` and `format = 1`. It also cannot be enabled with `flux_wrt`, `heat_ratio_wrt`, `pres_inf_wrt`, `c_wrt`, `omega_wrt`, `ib`, `schlieren_wrt`, `qm_wrt`, or 'liutex_wrt'. ### 8. Acoustic Source {#acoustic-source} diff --git a/examples/2D_lagrange_bubblescreen/case.py b/examples/2D_lagrange_bubblescreen/case.py index fb1dd1cf81..18fba87222 100644 --- a/examples/2D_lagrange_bubblescreen/case.py +++ b/examples/2D_lagrange_bubblescreen/case.py @@ -111,6 +111,7 @@ "precision": 2, "prim_vars_wrt": "T", "parallel_io": "T", + "lag_db_wrt": "T", # Patch 1: Water (left) "patch_icpp(1)%geometry": 3, "patch_icpp(1)%x_centroid": 0.0, diff --git a/examples/3D_lagrange_bubblescreen/case.py b/examples/3D_lagrange_bubblescreen/case.py index 1e0d81b81b..40a85eb022 100644 --- a/examples/3D_lagrange_bubblescreen/case.py +++ b/examples/3D_lagrange_bubblescreen/case.py @@ -120,6 +120,7 @@ "precision": 2, "prim_vars_wrt": "T", "parallel_io": "T", + "lag_db_wrt": "T", # Patch 1: Water (left) "patch_icpp(1)%geometry": 9, "patch_icpp(1)%x_centroid": 0.0, diff --git a/src/common/include/3dHardcodedIC.fpp b/src/common/include/3dHardcodedIC.fpp index 39fb581c0d..d23db26630 100644 --- a/src/common/include/3dHardcodedIC.fpp +++ b/src/common/include/3dHardcodedIC.fpp @@ -7,6 +7,67 @@ ! Case 302 - IGR Jet real(wp) :: r, ux_th, ux_am, p_th, p_am, rho_th, rho_am, y_th, z_th, r_th, eps_smooth + real(wp), allocatable, dimension(:, :) :: ih + integer :: pos, start, end + logical :: file_exist + character(len=10000) :: line + character(len=25) :: value + + if (patch_icpp(patch_id)%hcid == 303) then + allocate(ih(0:m_glb, 0:p_glb)) + + if (interface_file == '.') then + call s_mpi_abort("Error: interface_file must be specified for hcid=303") + else + inquire (file=trim(interface_file), exist=file_exist) + if (file_exist) then + open(unit=10, file=trim(interface_file), status="old", action="read") + do i = 0, m_glb + read(10, '(A)') line ! Read a full line as a string + start = 1 + + do j = 0, p_glb + end = index(line(start:), ',') ! Find the next comma + if (end == 0) then + value = trim(adjustl(line(start:))) ! Last value in the line + else + value = trim(adjustl(line(start:start+end-2))) ! Extract substring + start = start + end ! Move to next value + end if + read(value, *) ih(i, j) ! Convert string to numeric value + if (.not. f_is_default(normMag)) ih(i, j) = ih(i, j) * normMag + if (.not. f_is_default(normFac)) ih(i, j) = ih(i, j) + normFac + end do + end do + close(10) + else + call s_mpi_abort("Error: interface_file specified for hcid=303 does not exist") + end if + end if + end if + + if (patch_icpp(patch_id)%hcid == 304) then + allocate(ih(0:n_glb, 0:0)) + if (interface_file == '.') then + call s_mpi_abort("Error: interface_file must be specified for hcid=304") + else + inquire (file=trim(interface_file), exist=file_exist) + if (file_exist) then + open(unit=10, file=trim(interface_file), status="old", action="read") + do i = 0, n_glb + read(10, '(A)') line ! Read a full line as a string + value = trim(line) + read(value, *) ih(i, 0) ! Convert string to numeric value + if (.not. f_is_default(normMag)) ih(i, 0) = ih(i, 0) * normMag + if (.not. f_is_default(normFac)) ih(i, 0) = ih(i, 0) + normFac + end do + close(10) + else + call s_mpi_abort("Error: interface_file specified for hcid=304 does not exist") + end if + end if + end if + eps = 1e-9_wp #:enddef @@ -86,6 +147,40 @@ q_prim_vf(E_idx)%sf(i, j, k) = p_th*f_cut_on(r - r_th, eps_smooth)*f_cut_on(x_cc(i), eps_smooth) + p_am + case (303) ! 3D Interface from file cartesian + + alph = 0.5_wp * (1 + (1._wp - 2._wp * eps) * & + tanh((ih(start_idx(1) + i,start_idx(3) + k) - y_cc(j))*(0.5_wp / dx))) + + q_prim_vf(advxb)%sf(i,j,k) = alph + q_prim_vf(advxe)%sf(i,j,k) = 1._wp - alph + + q_prim_vf(contxb)%sf(i,j,k) = q_prim_vf(advxb)%sf(i,j,k) * 1._wp + q_prim_vf(contxe)%sf(i,j,k) = q_prim_vf(advxe)%sf(i,j,k) * (1._wp / 950._wp) + + q_prim_vf(E_idx)%sf(i,j,k) = p0 + & + (q_prim_vf(contxb)%sf(i,j,k) + q_prim_vf(contxe)%sf(i,j,k)) * g0 * & + (ih(start_idx(1) + i, start_idx(3) + k) - y_cc(j)) + + if (surface_tension) q_prim_vf(c_idx)%sf(i,j,k) = alph + + case (304) ! 3D Interface from file axisymmetric + + alph = 0.5_wp * (1 + (1._wp - 2._wp * eps) * & + tanh((ih(start_idx(2) + j,0) - x_cc(i))*(0.01_wp / dx))) + + q_prim_vf(advxb)%sf(i,j,k) = alph + q_prim_vf(advxe)%sf(i,j,k) = 1._wp - alph + + q_prim_vf(contxb)%sf(i,j,k) = q_prim_vf(advxb)%sf(i,j,k) * 1._wp + q_prim_vf(contxe)%sf(i,j,k) = q_prim_vf(advxe)%sf(i,j,k) * (1._wp / 950._wp) + + q_prim_vf(E_idx)%sf(i,j,k) = p0 + & + (q_prim_vf(contxb)%sf(i,j,k) + q_prim_vf(contxe)%sf(i,j,k)) * g0 * & + (ih(start_idx(1) + i, start_idx(3) + k) - y_cc(j)) + + if (surface_tension) q_prim_vf(c_idx)%sf(i,j,k) = alph + case (370) ! This hardcoded case extrudes a 2D profile to initialize a 3D simulation domain @: HardcodedReadValues() diff --git a/src/common/m_boundary_common.fpp b/src/common/m_boundary_common.fpp index f89278a86d..ee90fc7a61 100644 --- a/src/common/m_boundary_common.fpp +++ b/src/common/m_boundary_common.fpp @@ -1942,6 +1942,21 @@ contains offset_x%beg = buff_size; offset_x%end = buff_size offset_y%beg = buff_size; offset_y%end = buff_size offset_z%beg = buff_size; offset_z%end = buff_size + +#ifdef MFC_MPI + ! Populate global domain boundaries with stretched grids + call s_mpi_allreduce_min(x_cb(-1), glb_bounds(1)%beg) + call s_mpi_allreduce_max(x_cb(m), glb_bounds(1)%end) + + if (n > 0) then + call s_mpi_allreduce_min(y_cb(-1), glb_bounds(2)%beg) + call s_mpi_allreduce_max(y_cb(n), glb_bounds(2)%end) + if (p > 0) then + call s_mpi_allreduce_min(z_cb(-1), glb_bounds(3)%beg) + call s_mpi_allreduce_max(z_cb(p), glb_bounds(3)%end) + end if + end if +#endif #endif #ifndef MFC_PRE_PROCESS diff --git a/src/common/m_constants.fpp b/src/common/m_constants.fpp index 1d2e53d206..9fcd9c7134 100644 --- a/src/common/m_constants.fpp +++ b/src/common/m_constants.fpp @@ -13,6 +13,7 @@ module m_constants real(wp), parameter :: small_alf = 1.e-11_wp !< Small alf tolerance real(wp), parameter :: pi = 3.141592653589793_wp !< Pi real(wp), parameter :: verysmall = 1.e-12_wp !< Very small number + real(wp), parameter :: Re_b_min = 1.e-6_wp !< minimum bubble reynolds number for drag coeff calc. integer, parameter :: num_stcls_min = 5 !< Minimum # of stencils integer, parameter :: path_len = 400 !< Maximum path length @@ -21,7 +22,7 @@ module m_constants integer, parameter :: fourier_rings = 5 !< Fourier filter ring limit integer, parameter :: num_fluids_max = 10 !< Maximum number of fluids in the simulation integer, parameter :: num_probes_max = 10 !< Maximum number of flow probes in the simulation - integer, parameter :: num_patches_max = 10 + integer, parameter :: num_patches_max = 20 integer, parameter :: num_bc_patches_max = 10 integer, parameter :: pathlen_max = 400 integer, parameter :: nnode = 4 !< Number of QBMM nodes @@ -61,6 +62,7 @@ module m_constants ! Lagrange bubbles constants integer, parameter :: mapCells = 3 !< Number of cells around the bubble where the smoothening function will have effect real(wp), parameter :: R_uni = 8314._wp !< Universal gas constant - J/kmol/K + integer, parameter :: lag_io_vars = 21 ! Number of variables per particle for MPI_IO ! Strang Splitting constants real(wp), parameter :: dflt_adap_dt_tol = 1.e-4_wp !< Default tolerance for adaptive step size diff --git a/src/common/m_derived_types.fpp b/src/common/m_derived_types.fpp index 62c5fc9ceb..b529e1e2e5 100644 --- a/src/common/m_derived_types.fpp +++ b/src/common/m_derived_types.fpp @@ -447,6 +447,12 @@ module m_derived_types logical :: write_bubbles !< Write files to track the bubble evolution each time step logical :: write_bubbles_stats !< Write the maximum and minimum radius of each bubble integer :: nBubs_glb !< Global number of bubbles + integer :: vel_model !< Particle velocity model + integer :: drag_model !< Particle drag model + logical :: pressure_force !< Include pressure force translational motion + logical :: gravity_force !< Include gravity force in translational motion + logical :: momentum_transfer_force !< Include momentum transfer from radial dynamics in translational motion + real(wp) :: c_d !< Drag coefficient real(wp) :: epsilonb !< Standard deviation scaling for the gaussian function real(wp) :: charwidth !< Domain virtual depth (z direction, for 2D simulations) real(wp) :: valmaxvoid !< Maximum void fraction permitted diff --git a/src/common/m_helper_basic.fpp b/src/common/m_helper_basic.fpp index 697e919077..61a8d92849 100644 --- a/src/common/m_helper_basic.fpp +++ b/src/common/m_helper_basic.fpp @@ -112,10 +112,11 @@ contains pure subroutine s_configure_coordinate_bounds(recon_type, weno_polyn, muscl_polyn, & igr_order, buff_size, idwint, idwbuff, & - viscous, bubbles_lagrange, m, n, p, num_dims, igr, ib) + viscous, bubbles_lagrange, m, n, p, num_dims, & + igr, ib, fd_number) integer, intent(in) :: recon_type, weno_polyn, muscl_polyn - integer, intent(in) :: m, n, p, num_dims, igr_order + integer, intent(in) :: m, n, p, num_dims, igr_order, fd_number integer, intent(inout) :: buff_size type(int_bounds_info), dimension(3), intent(inout) :: idwint, idwbuff logical, intent(in) :: viscous, bubbles_lagrange @@ -140,7 +141,7 @@ contains ! Correction for smearing function in the lagrangian subgrid bubble model if (bubbles_lagrange) then - buff_size = max(buff_size, 6) + buff_size = max(buff_size + fd_number, 6 + fd_number) end if if (ib) then diff --git a/src/common/m_mpi_common.fpp b/src/common/m_mpi_common.fpp index 07463137f2..fd9310abc3 100644 --- a/src/common/m_mpi_common.fpp +++ b/src/common/m_mpi_common.fpp @@ -434,21 +434,40 @@ contains mpi_p, MPI_MIN, 0, & MPI_COMM_WORLD, ierr) end if - #else - icfl_max_glb = icfl_max_loc if (viscous) then vcfl_max_glb = vcfl_max_loc Rc_min_glb = Rc_min_loc end if - #endif #endif end subroutine s_mpi_reduce_stability_criteria_extrema + !> The following subroutine takes the inputted variable and + !! determines its sum on the entire computational domain. + !! @param var_loc holds the local value to be reduced among + !! all the processors in communicator. On output, the variable holds + !! the sum, reduced amongst all of the local values. + subroutine s_mpi_reduce_int_sum(var_loc, sum) + + integer, intent(inout) :: var_loc + integer, intent(inout) :: sum + +#ifdef MFC_MPI + integer :: ierr !< Generic flag used to identify and report MPI errors + + ! Performing reduction procedure and eventually storing its result + ! into the variable that was initially inputted into the subroutine + call MPI_REDUCE(var_loc, sum, 1, MPI_INTEGER, & + MPI_SUM, 0, MPI_COMM_WORLD, ierr) + +#endif + + end subroutine s_mpi_reduce_int_sum + !> The following subroutine takes the input local variable !! from all processors and reduces to the sum of all !! values. The reduced variable is recorded back onto the @@ -1148,9 +1167,17 @@ contains integer :: recon_order !< !! WENO or MUSCL reconstruction order - integer :: i, j !< Generic loop iterators + integer :: i, j, k !< Generic loop iterators integer :: ierr !< Generic flag used to identify and report MPI errors + ! temp array to store neighbor rank coordinates + integer, dimension(1:num_dims) :: neighbor_coords + + ! Zeroing out communication needs for moving EL bubbles/particles + nidx(1)%beg = 0; nidx(1)%end = 0 + nidx(2)%beg = 0; nidx(2)%end = 0 + nidx(3)%beg = 0; nidx(3)%end = 0 + if (recon_type == WENO_TYPE) then recon_order = weno_order else @@ -1327,6 +1354,7 @@ contains call MPI_CART_RANK(MPI_COMM_CART, proc_coords, & bc_z%beg, ierr) proc_coords(3) = proc_coords(3) + 1 + nidx(3)%beg = -1 end if ! Boundary condition at the end @@ -1335,6 +1363,7 @@ contains call MPI_CART_RANK(MPI_COMM_CART, proc_coords, & bc_z%end, ierr) proc_coords(3) = proc_coords(3) - 1 + nidx(3)%end = 1 end if #ifdef MFC_POST_PROCESS @@ -1466,6 +1495,7 @@ contains call MPI_CART_RANK(MPI_COMM_CART, proc_coords, & bc_y%beg, ierr) proc_coords(2) = proc_coords(2) + 1 + nidx(2)%beg = -1 end if ! Boundary condition at the end @@ -1474,6 +1504,7 @@ contains call MPI_CART_RANK(MPI_COMM_CART, proc_coords, & bc_y%end, ierr) proc_coords(2) = proc_coords(2) - 1 + nidx(2)%end = 1 end if #ifdef MFC_POST_PROCESS @@ -1559,6 +1590,7 @@ contains proc_coords(1) = proc_coords(1) - 1 call MPI_CART_RANK(MPI_COMM_CART, proc_coords, bc_x%beg, ierr) proc_coords(1) = proc_coords(1) + 1 + nidx(1)%beg = -1 end if ! Boundary condition at the end @@ -1566,18 +1598,19 @@ contains proc_coords(1) = proc_coords(1) + 1 call MPI_CART_RANK(MPI_COMM_CART, proc_coords, bc_x%end, ierr) proc_coords(1) = proc_coords(1) - 1 + nidx(1)%end = 1 end if #ifdef MFC_POST_PROCESS ! Ghost zone at the beginning - if (proc_coords(1) > 0 .and. format == 1 .and. n > 0) then + if (proc_coords(1) > 0 .and. format == 1) then offset_x%beg = 2 else offset_x%beg = 0 end if ! Ghost zone at the end - if (proc_coords(1) < num_procs_x - 1 .and. format == 1 .and. n > 0) then + if (proc_coords(1) < num_procs_x - 1 .and. format == 1) then offset_x%end = 2 else offset_x%end = 0 @@ -1611,6 +1644,24 @@ contains end if #endif end if + + @:ALLOCATE(neighbor_ranks(nidx(1)%beg:nidx(1)%end, & + nidx(2)%beg:nidx(2)%end, & + nidx(3)%beg:nidx(3)%end)) + + do k = nidx(3)%beg, nidx(3)%end + do j = nidx(2)%beg, nidx(2)%end + do i = nidx(1)%beg, nidx(1)%end + if (abs(i) + abs(j) + abs(k) > 0) then + neighbor_coords(1) = proc_coords(1) + i + if (num_dims > 1) neighbor_coords(2) = proc_coords(2) + j + if (num_dims > 2) neighbor_coords(3) = proc_coords(3) + k + call MPI_CART_RANK(MPI_COMM_CART, neighbor_coords, & + neighbor_ranks(i, j, k), ierr) + end if + end do + end do + end do #endif end subroutine s_mpi_decompose_computational_domain diff --git a/src/post_process/m_data_output.fpp b/src/post_process/m_data_output.fpp index 4c8867225e..ea901fa438 100644 --- a/src/post_process/m_data_output.fpp +++ b/src/post_process/m_data_output.fpp @@ -30,7 +30,8 @@ module m_data_output s_open_energy_data_file, & s_write_grid_to_formatted_database_file, & s_write_variable_to_formatted_database_file, & - s_write_lag_bubbles_results, & + s_write_lag_bubbles_results_to_text, & + s_write_lag_bubbles_to_formatted_database_file, & s_write_intf_data_file, & s_write_energy_data_file, & s_close_formatted_database_file, & @@ -155,7 +156,7 @@ contains ! the offsets and the one bookkeeping the number of cell-boundaries ! in each active coordinate direction. Note that all these variables ! are only needed by the Silo-HDF5 format for multidimensional data. - if (format == 1 .and. n > 0) then + if (format == 1) then allocate (data_extents(1:2, 0:num_procs - 1)) @@ -164,11 +165,16 @@ contains allocate (lo_offset(1:3)) allocate (hi_offset(1:3)) allocate (dims(1:3)) - else + elseif (n > 0) then allocate (spatial_extents(1:4, 0:num_procs - 1)) allocate (lo_offset(1:2)) allocate (hi_offset(1:2)) allocate (dims(1:2)) + else + allocate (spatial_extents(1:2, 0:num_procs - 1)) + allocate (lo_offset(1:1)) + allocate (hi_offset(1:1)) + allocate (dims(1:1)) end if end if @@ -180,7 +186,7 @@ contains ! With the same, latter, requirements, the variables bookkeeping the ! number of cell-boundaries in each active coordinate direction are ! also set here. - if (format == 1 .and. n > 0) then + if (format == 1) then if (p > 0) then if (grid_geometry == 3) then lo_offset(:) = (/offset_y%beg, offset_z%beg, offset_x%beg/) @@ -199,12 +205,16 @@ contains n + offset_y%beg + offset_y%end + 2, & p + offset_z%beg + offset_z%end + 2/) end if - else + elseif (n > 0) then lo_offset(:) = (/offset_x%beg, offset_y%beg/) hi_offset(:) = (/offset_x%end, offset_y%end/) dims(:) = (/m + offset_x%beg + offset_x%end + 2, & n + offset_y%beg + offset_y%end + 2/) + else + lo_offset(:) = (/offset_x%beg/) + hi_offset(:) = (/offset_x%end/) + dims(:) = (/m + offset_x%beg + offset_x%end + 2/) end if end if @@ -276,12 +286,14 @@ contains end if if (bubbles_lagrange) then !Lagrangian solver - dbdir = trim(case_dir)//'/lag_bubbles_post_process' - file_loc = trim(dbdir)//'/.' - call my_inquire(file_loc, dir_check) + if (lag_txt_wrt) then + dbdir = trim(case_dir)//'/lag_bubbles_post_process' + file_loc = trim(dbdir)//'/.' + call my_inquire(file_loc, dir_check) - if (dir_check .neqv. .true.) then - call s_create_directory(trim(dbdir)) + if (dir_check .neqv. .true.) then + call s_create_directory(trim(dbdir)) + end if end if end if @@ -656,7 +668,7 @@ contains ! Silo-HDF5 Database Format - if (format == 1 .and. n > 0) then + if (format == 1) then ! For multidimensional data sets, the spatial extents of all of ! the grid(s) handled by the local processor(s) are recorded so @@ -664,7 +676,6 @@ contains ! database master file. if (num_procs > 1) then call s_mpi_gather_spatial_extents(spatial_extents) - elseif (p > 0) then if (grid_geometry == 3) then spatial_extents(:, 0) = (/minval(y_cb), minval(z_cb), & @@ -675,11 +686,11 @@ contains minval(z_cb), maxval(x_cb), & maxval(y_cb), maxval(z_cb)/) end if - - else + elseif (n > 0) then spatial_extents(:, 0) = (/minval(x_cb), minval(y_cb), & maxval(x_cb), maxval(y_cb)/) - + else + spatial_extents(:, 0) = (/minval(x_cb), maxval(x_cb)/) end if ! Next, the root processor proceeds to record all of the spatial @@ -712,48 +723,45 @@ contains ! with its offsets that indicate the presence and size of ghost ! zone layer(s), are put in the formatted database slave file. - if (precision == 1) then - if (p > 0) then - z_cb_s(:) = real(z_cb(:), sp) + if (p > 0) then + err = DBMKOPTLIST(2, optlist) + err = DBADDIOPT(optlist, DBOPT_LO_OFFSET, lo_offset) + err = DBADDIOPT(optlist, DBOPT_HI_OFFSET, hi_offset) + if (grid_geometry == 3) then + err = DBPUTQM(dbfile, 'rectilinear_grid', 16, & + 'x', 1, 'y', 1, 'z', 1, & + y_cb, z_cb, x_cb, dims, 3, & + DB_DOUBLE, DB_COLLINEAR, & + optlist, ierr) + else + err = DBPUTQM(dbfile, 'rectilinear_grid', 16, & + 'x', 1, 'y', 1, 'z', 1, & + x_cb, y_cb, z_cb, dims, 3, & + DB_DOUBLE, DB_COLLINEAR, & + optlist, ierr) end if - x_cb_s(:) = real(x_cb(:), sp) - y_cb_s(:) = real(y_cb(:), sp) + err = DBFREEOPTLIST(optlist) + elseif (n > 0) then + err = DBMKOPTLIST(2, optlist) + err = DBADDIOPT(optlist, DBOPT_LO_OFFSET, lo_offset) + err = DBADDIOPT(optlist, DBOPT_HI_OFFSET, hi_offset) + err = DBPUTQM(dbfile, 'rectilinear_grid', 16, & + 'x', 1, 'y', 1, 'z', 1, & + x_cb, y_cb, DB_F77NULL, dims, 2, & + DB_DOUBLE, DB_COLLINEAR, & + optlist, ierr) + err = DBFREEOPTLIST(optlist) + else + err = DBMKOPTLIST(2, optlist) + err = DBADDIOPT(optlist, DBOPT_LO_OFFSET, lo_offset) + err = DBADDIOPT(optlist, DBOPT_HI_OFFSET, hi_offset) + err = DBPUTQM(dbfile, 'rectilinear_grid', 16, & + 'x', 1, 'y', 1, 'z', 1, & + x_cb, DB_F77NULL, DB_F77NULL, dims, 1, & + DB_DOUBLE, DB_COLLINEAR, & + optlist, ierr) + err = DBFREEOPTLIST(optlist) end if - - #:for PRECISION, SFX, DBT in [(1,'_s','DB_FLOAT'),(2,'',"DB_DOUBLE")] - if (precision == ${PRECISION}$) then - if (p > 0) then - err = DBMKOPTLIST(2, optlist) - err = DBADDIOPT(optlist, DBOPT_LO_OFFSET, lo_offset) - err = DBADDIOPT(optlist, DBOPT_HI_OFFSET, hi_offset) - if (grid_geometry == 3) then - err = DBPUTQM(dbfile, 'rectilinear_grid', 16, & - 'x', 1, 'y', 1, 'z', 1, & - y_cb${SFX}$, z_cb${SFX}$, x_cb${SFX}$, dims, 3, & - ${DBT}$, DB_COLLINEAR, & - optlist, ierr) - else - err = DBPUTQM(dbfile, 'rectilinear_grid', 16, & - 'x', 1, 'y', 1, 'z', 1, & - x_cb${SFX}$, y_cb${SFX}$, z_cb${SFX}$, dims, 3, & - ${DBT}$, DB_COLLINEAR, & - optlist, ierr) - end if - err = DBFREEOPTLIST(optlist) - else - err = DBMKOPTLIST(2, optlist) - err = DBADDIOPT(optlist, DBOPT_LO_OFFSET, lo_offset) - err = DBADDIOPT(optlist, DBOPT_HI_OFFSET, hi_offset) - err = DBPUTQM(dbfile, 'rectilinear_grid', 16, & - 'x', 1, 'y', 1, 'z', 1, & - x_cb${SFX}$, y_cb${SFX}$, DB_F77NULL, dims, 2, & - ${DBT}$, DB_COLLINEAR, & - optlist, ierr) - err = DBFREEOPTLIST(optlist) - end if - end if - #:endfor - ! END: Silo-HDF5 Database Format ! Binary Database Format @@ -870,144 +878,46 @@ contains if (format == 1) then - ! In 1D, a curve object, featuring the local processor grid and - ! flow variable data, is written to the formatted database slave - ! file. The root process, on the other hand, will also take care - ! of gathering the entire grid and associated flow variable data - ! and write it to the formatted database master file. - if (n == 0) then - - if (precision == 1 .and. wp == dp) then - x_cc_s(:) = real(x_cc(:), sp) - q_sf_s(:, :, :) = real(q_sf(:, :, :), sp) - elseif (precision == 1 .and. wp == sp) then - x_cc_s(:) = x_cc(:) - q_sf_s(:, :, :) = q_sf(:, :, :) - end if - - ! Writing the curve object associated with the local process - ! to the formatted database slave file - #:for PRECISION, SFX, DBT in [(1,'_s','DB_FLOAT'),(2,'',"DB_DOUBLE")] - if (precision == ${PRECISION}$) then - err = DBPUTCURVE(dbfile, trim(varname), len_trim(varname), & - x_cc${SFX}$ (0:m), q_sf${SFX}$, ${DBT}$, m + 1, & - DB_F77NULL, ierr) - end if - #:endfor - - ! Assembling the local grid and flow variable data for the - ! entire computational domain on to the root process - - if (num_procs > 1) then - call s_mpi_defragment_1d_grid_variable() - call s_mpi_defragment_1d_flow_variable(q_sf, q_root_sf) - - if (precision == 1) then - x_root_cc_s(:) = real(x_root_cc(:), sp) - q_root_sf_s(:, :, :) = real(q_root_sf(:, :, :), sp) - end if - else - if (precision == 1) then - x_root_cc_s(:) = real(x_cc(:), sp) - q_root_sf_s(:, :, :) = real(q_sf(:, :, :), sp) - else - x_root_cc(:) = x_cc(:) - q_root_sf(:, :, :) = q_sf(:, :, :) - end if - end if - - ! Writing the curve object associated with the root process - ! to the formatted database master file - if (proc_rank == 0) then - #:for PRECISION, SFX, DBT in [(1,'_s','DB_FLOAT'),(2,'',"DB_DOUBLE")] - if (precision == ${PRECISION}$) then - err = DBPUTCURVE(dbroot, trim(varname), & - len_trim(varname), & - x_root_cc${SFX}$, q_root_sf${SFX}$, & - ${DBT}$, m_root + 1, & - DB_F77NULL, ierr) - end if - #:endfor - end if - - return - - ! In multidimensions, the local process(es) take care of writing - ! the flow variable data they are in charge of to the formatted - ! database slave file. The root processor, additionally, is also - ! responsible in gathering the flow variable extents of each of - ! the local processor(s) and writing them to formatted database - ! master file. + ! Determining the extents of the flow variable on each local + ! process and gathering all this information on root process + if (num_procs > 1) then + call s_mpi_gather_data_extents(q_sf, data_extents) else + data_extents(:, 0) = (/minval(q_sf), maxval(q_sf)/) + end if - ! Determining the extents of the flow variable on each local - ! process and gathering all this information on root process - if (num_procs > 1) then - call s_mpi_gather_data_extents(q_sf, data_extents) - else - data_extents(:, 0) = (/minval(q_sf), maxval(q_sf)/) - end if - - ! Next, the root process proceeds to write the gathered flow - ! variable data extents to formatted database master file. - if (proc_rank == 0) then + ! Next, the root process proceeds to write the gathered flow + ! variable data extents to formatted database master file. + if (proc_rank == 0) then - do i = 1, num_procs - write (varnames(i), '(A,I0,A,I0,A)') '../p', i - 1, & - '/', t_step, '.silo:'//trim(varname) - end do + do i = 1, num_procs + write (varnames(i), '(A,I0,A,I0,A)') '../p', i - 1, & + '/', t_step, '.silo:'//trim(varname) + end do - vartypes = DB_QUADVAR + vartypes = DB_QUADVAR - err = DBSET2DSTRLEN(len(varnames(1))) - err = DBMKOPTLIST(2, optlist) - err = DBADDIOPT(optlist, DBOPT_EXTENTS_SIZE, 2) - err = DBADDDOPT(optlist, DBOPT_EXTENTS, data_extents) - err = DBPUTMVAR(dbroot, trim(varname), & - len_trim(varname), num_procs, & - varnames, len_trim(varnames), & - vartypes, optlist, ierr) - err = DBFREEOPTLIST(optlist) + err = DBSET2DSTRLEN(len(varnames(1))) + err = DBMKOPTLIST(2, optlist) + err = DBADDIOPT(optlist, DBOPT_EXTENTS_SIZE, 2) + err = DBADDDOPT(optlist, DBOPT_EXTENTS, data_extents) + err = DBPUTMVAR(dbroot, trim(varname), & + len_trim(varname), num_procs, & + varnames, len_trim(varnames), & + vartypes, optlist, ierr) + err = DBFREEOPTLIST(optlist) - end if + end if - ! Finally, each of the local processor(s) proceeds to write - ! the flow variable data that it is responsible for to the - ! formatted database slave file. - if (wp == dp) then - if (precision == 1) then - do i = -offset_x%beg, m + offset_x%end - do j = -offset_y%beg, n + offset_y%end - do k = -offset_z%beg, p + offset_z%end - q_sf_s(i, j, k) = real(q_sf(i, j, k), sp) - end do - end do - end do - if (grid_geometry == 3) then - do i = -offset_x%beg, m + offset_x%end - do j = -offset_y%beg, n + offset_y%end - do k = -offset_z%beg, p + offset_z%end - cyl_q_sf_s(j, k, i) = q_sf_s(i, j, k) - end do - end do - end do - end if - else - if (grid_geometry == 3) then - do i = -offset_x%beg, m + offset_x%end - do j = -offset_y%beg, n + offset_y%end - do k = -offset_z%beg, p + offset_z%end - cyl_q_sf(j, k, i) = q_sf(i, j, k) - end do - end do - end do - end if - end if - elseif (wp == sp) then + ! Finally, each of the local processor(s) proceeds to write + ! the flow variable data that it is responsible for to the + ! formatted database slave file. + if (wp == dp) then + if (precision == 1) then do i = -offset_x%beg, m + offset_x%end do j = -offset_y%beg, n + offset_y%end do k = -offset_z%beg, p + offset_z%end - q_sf_s(i, j, k) = q_sf(i, j, k) + q_sf_s(i, j, k) = real(q_sf(i, j, k), sp) end do end do end do @@ -1020,38 +930,72 @@ contains end do end do end if + else + if (grid_geometry == 3) then + do i = -offset_x%beg, m + offset_x%end + do j = -offset_y%beg, n + offset_y%end + do k = -offset_z%beg, p + offset_z%end + cyl_q_sf(j, k, i) = q_sf(i, j, k) + end do + end do + end do + end if + end if + elseif (wp == sp) then + do i = -offset_x%beg, m + offset_x%end + do j = -offset_y%beg, n + offset_y%end + do k = -offset_z%beg, p + offset_z%end + q_sf_s(i, j, k) = q_sf(i, j, k) + end do + end do + end do + if (grid_geometry == 3) then + do i = -offset_x%beg, m + offset_x%end + do j = -offset_y%beg, n + offset_y%end + do k = -offset_z%beg, p + offset_z%end + cyl_q_sf_s(j, k, i) = q_sf_s(i, j, k) + end do + end do + end do end if + end if - #:for PRECISION, SFX, DBT in [(1,'_s','DB_FLOAT'),(2,'',"DB_DOUBLE")] - if (precision == ${PRECISION}$) then - if (p > 0) then - if (grid_geometry == 3) then - err = DBPUTQV1(dbfile, trim(varname), & - len_trim(varname), & - 'rectilinear_grid', 16, & - cyl_q_sf${SFX}$, dims - 1, 3, DB_F77NULL, & - 0, ${DBT}$, DB_ZONECENT, & - DB_F77NULL, ierr) - else - err = DBPUTQV1(dbfile, trim(varname), & - len_trim(varname), & - 'rectilinear_grid', 16, & - q_sf${SFX}$, dims - 1, 3, DB_F77NULL, & - 0, ${DBT}$, DB_ZONECENT, & - DB_F77NULL, ierr) - end if + #:for PRECISION, SFX, DBT in [(1,'_s','DB_FLOAT'),(2,'',"DB_DOUBLE")] + if (precision == ${PRECISION}$) then + if (p > 0) then + if (grid_geometry == 3) then + err = DBPUTQV1(dbfile, trim(varname), & + len_trim(varname), & + 'rectilinear_grid', 16, & + cyl_q_sf${SFX}$, dims - 1, 3, DB_F77NULL, & + 0, ${DBT}$, DB_ZONECENT, & + DB_F77NULL, ierr) else err = DBPUTQV1(dbfile, trim(varname), & len_trim(varname), & 'rectilinear_grid', 16, & - q_sf${SFX}$, dims - 1, 2, DB_F77NULL, & + q_sf${SFX}$, dims - 1, 3, DB_F77NULL, & 0, ${DBT}$, DB_ZONECENT, & DB_F77NULL, ierr) end if - end if - #:endfor + elseif (n > 0) then + err = DBPUTQV1(dbfile, trim(varname), & + len_trim(varname), & + 'rectilinear_grid', 16, & + q_sf${SFX}$, dims - 1, 2, DB_F77NULL, & + 0, ${DBT}$, DB_ZONECENT, & + DB_F77NULL, ierr) + else + err = DBPUTQV1(dbfile, trim(varname), & + len_trim(varname), & + 'rectilinear_grid', 16, & + q_sf${SFX}$, dims - 1, 1, DB_F77NULL, & + 0, ${DBT}$, DB_ZONECENT, & + DB_F77NULL, ierr) - end if + end if + end if + #:endfor ! END: Silo-HDF5 Database Format @@ -1094,7 +1038,7 @@ contains !> Subroutine that writes the post processed results in the folder 'lag_bubbles_data' !! @param t_step Current time step - impure subroutine s_write_lag_bubbles_results(t_step) + impure subroutine s_write_lag_bubbles_results_to_text(t_step) integer, intent(in) :: t_step @@ -1112,32 +1056,66 @@ contains logical :: lg_bub_file, file_exist integer, dimension(2) :: gsizes, lsizes, start_idx_part - integer :: ifile, tot_data + integer :: ifile integer :: ierr !< Generic flag used to identify and report MPI errors + real(wp) :: file_time, file_dt + integer :: file_num_procs, file_tot_part, tot_part integer :: i - write (file_loc, '(A,I0,A)') 'lag_bubbles_mpi_io_', t_step, '.dat' + integer, dimension(:), allocatable :: proc_bubble_counts + real(wp), dimension(1:1, 1:lag_io_vars) :: lag_io_null + lag_io_null = 0._wp + + ! Construct file path + write (file_loc, '(A,I0,A)') 'lag_bubbles_', t_step, '.dat' file_loc = trim(case_dir)//'/restart_data'//trim(mpiiofs)//trim(file_loc) + + ! Check if file exists inquire (FILE=trim(file_loc), EXIST=file_exist) + if (.not. file_exist) then + call s_mpi_abort('Restart file '//trim(file_loc)//' does not exist!') + end if - if (file_exist) then - if (proc_rank == 0) then - open (9, FILE=trim(file_loc), FORM='unformatted', STATUS='unknown') - read (9) tot_data, time_real - close (9) - end if - else - print '(A)', trim(file_loc)//' is missing. Exiting.' - call s_mpi_abort + if (.not. parallel_io) return + + if (proc_rank == 0) then + call MPI_FILE_OPEN(MPI_COMM_SELF, file_loc, MPI_MODE_RDONLY, & + mpi_info_int, ifile, ierr) + + call MPI_FILE_READ(ifile, file_tot_part, 1, MPI_INTEGER, status, ierr) + call MPI_FILE_READ(ifile, file_time, 1, mpi_p, status, ierr) + call MPI_FILE_READ(ifile, file_dt, 1, mpi_p, status, ierr) + call MPI_FILE_READ(ifile, file_num_procs, 1, MPI_INTEGER, status, ierr) + + call MPI_FILE_CLOSE(ifile, ierr) + end if + + call MPI_BCAST(file_tot_part, 1, MPI_INTEGER, 0, MPI_COMM_WORLD, ierr) + call MPI_BCAST(file_time, 1, mpi_p, 0, MPI_COMM_WORLD, ierr) + call MPI_BCAST(file_dt, 1, mpi_p, 0, MPI_COMM_WORLD, ierr) + call MPI_BCAST(file_num_procs, 1, MPI_INTEGER, 0, MPI_COMM_WORLD, ierr) + + allocate (proc_bubble_counts(file_num_procs)) + + if (proc_rank == 0) then + call MPI_FILE_OPEN(MPI_COMM_SELF, file_loc, MPI_MODE_RDONLY, & + mpi_info_int, ifile, ierr) + + ! Skip to processor counts position + disp = int(sizeof(file_tot_part) + 2*sizeof(file_time) + sizeof(file_num_procs), & + MPI_OFFSET_KIND) + call MPI_FILE_SEEK(ifile, disp, MPI_SEEK_SET, ierr) + call MPI_FILE_READ(ifile, proc_bubble_counts, file_num_procs, MPI_INTEGER, status, ierr) + + call MPI_FILE_CLOSE(ifile, ierr) end if - call MPI_BCAST(tot_data, 1, MPI_integer, 0, MPI_COMM_WORLD, ierr) - call MPI_BCAST(time_real, 1, mpi_p, 0, MPI_COMM_WORLD, ierr) + call MPI_BCAST(proc_bubble_counts, file_num_procs, MPI_INTEGER, 0, MPI_COMM_WORLD, ierr) - gsizes(1) = tot_data - gsizes(2) = 21 - lsizes(1) = tot_data - lsizes(2) = 21 + gsizes(1) = file_tot_part + gsizes(2) = lag_io_vars + lsizes(1) = file_tot_part + lsizes(2) = lag_io_vars start_idx_part(1) = 0 start_idx_part(2) = 0 @@ -1145,59 +1123,372 @@ contains MPI_ORDER_FORTRAN, mpi_p, view, ierr) call MPI_TYPE_COMMIT(view, ierr) + call MPI_FILE_OPEN(MPI_COMM_WORLD, file_loc, MPI_MODE_RDONLY, & + mpi_info_int, ifile, ierr) + + disp = int(sizeof(file_tot_part) + 2*sizeof(file_time) + sizeof(file_num_procs) + & + file_num_procs*sizeof(proc_bubble_counts(1)), MPI_OFFSET_KIND) + call MPI_FILE_SET_VIEW(ifile, disp, mpi_p, view, & + 'native', mpi_info_null, ierr) + + allocate (MPI_IO_DATA_lg_bubbles(file_tot_part, 1:lag_io_vars)) + + call MPI_FILE_READ_ALL(ifile, MPI_IO_DATA_lg_bubbles, lag_io_vars*file_tot_part, & + mpi_p, status, ierr) + + write (file_loc, '(A,I0,A)') 'lag_bubbles_post_process_', t_step, '.dat' + file_loc = trim(case_dir)//'/lag_bubbles_post_process/'//trim(file_loc) + + if (proc_rank == 0) then + open (unit=29, file=file_loc, form='formatted', position='rewind') + + if (lag_header) then + write (29, '(A)', advance='no') + if (lag_id_wrt) write (29, '(A8)', advance='no') 'id, ' + if (lag_pos_wrt) write (29, '(3(A17))', advance='no') 'px, ', 'py, ', 'pz, ' + if (lag_pos_prev_wrt) write (29, '(3(A17))', advance='no') 'pvx, ', 'pvy, ', 'pvz, ' + if (lag_vel_wrt) write (29, '(3(A17))', advance='no') 'vx, ', 'vy, ', 'vz, ' + if (lag_rad_wrt) write (29, '(A17)', advance='no') 'radius, ' + if (lag_rvel_wrt) write (29, '(A17)', advance='no') 'rvel, ' + if (lag_r0_wrt) write (29, '(A17)', advance='no') 'r0, ' + if (lag_rmax_wrt) write (29, '(A17)', advance='no') 'rmax, ' + if (lag_rmin_wrt) write (29, '(A17)', advance='no') 'rmin, ' + if (lag_dphidt_wrt) write (29, '(A17)', advance='no') 'dphidt, ' + if (lag_pres_wrt) write (29, '(A17)', advance='no') 'pressure, ' + if (lag_mv_wrt) write (29, '(A17)', advance='no') 'mv, ' + if (lag_mg_wrt) write (29, '(A17)', advance='no') 'mg, ' + if (lag_betaT_wrt) write (29, '(A17)', advance='no') 'betaT, ' + if (lag_betaC_wrt) write (29, '(A17)', advance='no') 'betaC, ' + write (29, '(A15)') 'time' + end if + + do i = 1, file_tot_part + id = int(MPI_IO_DATA_lg_bubbles(i, 1)) + inputvals(1:20) = MPI_IO_DATA_lg_bubbles(i, 2:21) + if (id > 0) then + write (29, '(100(A))', advance='no') '' + if (lag_id_wrt) write (29, '(I6, A)', advance='no') id, ', ' + if (lag_pos_wrt) write (29, '(3(E15.7, A))', advance='no') inputvals(1), ', ', inputvals(2), ', ', inputvals(3), ', ' + if (lag_pos_prev_wrt) write (29, '(3(E15.7, A))', advance='no') inputvals(4), ', ', inputvals(5), ', ', inputvals(6), ', ' + if (lag_vel_wrt) write (29, '(3(E15.7, A))', advance='no') inputvals(7), ', ', inputvals(8), ', ', inputvals(9), ', ' + if (lag_rad_wrt) write (29, '(E15.7, A)', advance='no') inputvals(10), ', ' + if (lag_rvel_wrt) write (29, '(E15.7, A)', advance='no') inputvals(11), ', ' + if (lag_r0_wrt) write (29, '(E15.7, A)', advance='no') inputvals(12), ', ' + if (lag_rmax_wrt) write (29, '(E15.7, A)', advance='no') inputvals(13), ', ' + if (lag_rmin_wrt) write (29, '(E15.7, A)', advance='no') inputvals(14), ', ' + if (lag_dphidt_wrt) write (29, '(E15.7, A)', advance='no') inputvals(15), ', ' + if (lag_pres_wrt) write (29, '(E15.7, A)', advance='no') inputvals(16), ', ' + if (lag_mv_wrt) write (29, '(E15.7, A)', advance='no') inputvals(17), ', ' + if (lag_mg_wrt) write (29, '(E15.7, A)', advance='no') inputvals(18), ', ' + if (lag_betaT_wrt) write (29, '(E15.7, A)', advance='no') inputvals(19), ', ' + if (lag_betaC_wrt) write (29, '(E15.7, A)', advance='no') inputvals(20), ', ' + write (29, '(E15.7)') time_real + end if + end do + close (29) + end if + + deallocate (MPI_IO_DATA_lg_bubbles) + + call s_mpi_barrier() + + call MPI_FILE_CLOSE(ifile, ierr) +#endif + + end subroutine s_write_lag_bubbles_results_to_text + + impure subroutine s_write_lag_bubbles_to_formatted_database_file(t_step) + + integer, intent(in) :: t_step + + character(len=len_trim(case_dir) + 3*name_len) :: file_loc + + integer :: id + +#ifdef MFC_MPI + real(wp), dimension(20) :: inputvals + real(wp) :: time_real + integer, dimension(MPI_STATUS_SIZE) :: status + integer(KIND=MPI_OFFSET_KIND) :: disp + integer :: view + + logical :: lg_bub_file, file_exist + + integer, dimension(2) :: gsizes, lsizes, start_idx_part + integer :: ifile, ierr, tot_data, valid_data, nBub + real(wp) :: file_time, file_dt + integer :: file_num_procs, file_tot_part + integer, dimension(:), allocatable :: proc_bubble_counts + real(wp), dimension(1:1, 1:lag_io_vars) :: dummy + character(LEN=4*name_len), dimension(num_procs) :: meshnames + integer, dimension(num_procs) :: meshtypes + real(wp) :: dummy_data + + integer :: i, j + + real(wp), dimension(:), allocatable :: bub_id + real(wp), dimension(:), allocatable :: px, py, pz, ppx, ppy, ppz, vx, vy, vz + real(wp), dimension(:), allocatable :: radius, rvel, rnot, rmax, rmin, dphidt + real(wp), dimension(:), allocatable :: pressure, mv, mg, betaT, betaC + + dummy = 0._wp + dummy_data = 0._wp + + ! Construct file path write (file_loc, '(A,I0,A)') 'lag_bubbles_', t_step, '.dat' file_loc = trim(case_dir)//'/restart_data'//trim(mpiiofs)//trim(file_loc) - inquire (FILE=trim(file_loc), EXIST=lg_bub_file) - if (lg_bub_file) then + ! Check if file exists + inquire (FILE=trim(file_loc), EXIST=file_exist) + if (.not. file_exist) then + call s_mpi_abort('Restart file '//trim(file_loc)//' does not exist!') + end if + + if (.not. parallel_io) return + + if (proc_rank == 0) then + call MPI_FILE_OPEN(MPI_COMM_SELF, file_loc, MPI_MODE_RDONLY, & + mpi_info_int, ifile, ierr) + + call MPI_FILE_READ(ifile, file_tot_part, 1, MPI_INTEGER, status, ierr) + call MPI_FILE_READ(ifile, file_time, 1, mpi_p, status, ierr) + call MPI_FILE_READ(ifile, file_dt, 1, mpi_p, status, ierr) + call MPI_FILE_READ(ifile, file_num_procs, 1, MPI_INTEGER, status, ierr) + + call MPI_FILE_CLOSE(ifile, ierr) + end if + + call MPI_BCAST(file_tot_part, 1, MPI_INTEGER, 0, MPI_COMM_WORLD, ierr) + call MPI_BCAST(file_time, 1, mpi_p, 0, MPI_COMM_WORLD, ierr) + call MPI_BCAST(file_dt, 1, mpi_p, 0, MPI_COMM_WORLD, ierr) + call MPI_BCAST(file_num_procs, 1, MPI_INTEGER, 0, MPI_COMM_WORLD, ierr) + + allocate (proc_bubble_counts(file_num_procs)) + + if (proc_rank == 0) then + call MPI_FILE_OPEN(MPI_COMM_SELF, file_loc, MPI_MODE_RDONLY, & + mpi_info_int, ifile, ierr) + + ! Skip to processor counts position + disp = int(sizeof(file_tot_part) + 2*sizeof(file_time) + sizeof(file_num_procs), & + MPI_OFFSET_KIND) + call MPI_FILE_SEEK(ifile, disp, MPI_SEEK_SET, ierr) + call MPI_FILE_READ(ifile, proc_bubble_counts, file_num_procs, MPI_INTEGER, status, ierr) + + call MPI_FILE_CLOSE(ifile, ierr) + end if + + call MPI_BCAST(proc_bubble_counts, file_num_procs, MPI_INTEGER, 0, MPI_COMM_WORLD, ierr) + + ! Set time variables from file + + nBub = proc_bubble_counts(proc_rank + 1) + + start_idx_part(1) = 0 + do i = 1, proc_rank + start_idx_part(1) = start_idx_part(1) + proc_bubble_counts(i) + end do + + start_idx_part(2) = 0 + lsizes(1) = nBub + lsizes(2) = lag_io_vars + + gsizes(1) = file_tot_part + gsizes(2) = lag_io_vars + + if (nBub > 0) then + + #:for VAR in ['bub_id', 'px', 'py', 'pz', 'ppx', 'ppy', 'ppz', 'vx', 'vy', 'vz', & + 'radius', 'rvel', 'rnot', 'rmax', 'rmin', 'dphidt', & + 'pressure', 'mv', 'mg', 'betaT', 'betaC'] + allocate (${VAR}$ (nBub)) + #:endfor + allocate (MPI_IO_DATA_lg_bubbles(nBub, 1:lag_io_vars)) + + call MPI_TYPE_CREATE_SUBARRAY(2, gsizes, lsizes, start_idx_part, & + MPI_ORDER_FORTRAN, mpi_p, view, ierr) + call MPI_TYPE_COMMIT(view, ierr) call MPI_FILE_OPEN(MPI_COMM_WORLD, file_loc, MPI_MODE_RDONLY, & mpi_info_int, ifile, ierr) - disp = 0._wp - call MPI_FILE_SET_VIEW(ifile, disp, mpi_p, view, & - 'native', mpi_info_null, ierr) + ! Skip extended header + disp = int(sizeof(file_tot_part) + 2*sizeof(file_time) + sizeof(file_num_procs) + & + file_num_procs*sizeof(proc_bubble_counts(1)), MPI_OFFSET_KIND) + call MPI_FILE_SET_VIEW(ifile, disp, mpi_p, view, 'native', mpi_info_int, ierr) - allocate (MPI_IO_DATA_lg_bubbles(tot_data, 1:21)) + call MPI_FILE_READ_ALL(ifile, MPI_IO_DATA_lg_bubbles, & + lag_io_vars*nBub, mpi_p, status, ierr) - call MPI_FILE_READ_ALL(ifile, MPI_IO_DATA_lg_bubbles, 21*tot_data, & - mpi_p, status, ierr) + call MPI_FILE_CLOSE(ifile, ierr) + call MPI_TYPE_FREE(view, ierr) - write (file_loc, '(A,I0,A)') 'lag_bubbles_post_process_', t_step, '.dat' - file_loc = trim(case_dir)//'/lag_bubbles_post_process/'//trim(file_loc) + ! Extract data from MPI_IO_DATA_lg_bubbles array + ! Adjust these indices based on your actual data layout + #:for VAR, IDX in [('bub_id', 1), ('px', 2), ('py',3), ('pz',4), ('ppx',5), ('ppy',6), ('ppz',7), & + ('vx',8), ('vy',9), ('vz',10), ('radius',11), ('rvel',12), & + ('rnot',13), ('rmax',14), ('rmin',15), ('dphidt',16), & + ('pressure',17), ('mv',18), ('mg',19), ('betaT',20), ('betaC',21)] + ${VAR}$ (:) = MPI_IO_DATA_lg_bubbles(:, ${IDX}$) + #:endfor + ! Next, the root processor proceeds to record all of the spatial + ! extents in the formatted database master file. In addition, it + ! also records a sub-domain connectivity map so that the entire + ! grid may be reassembled by looking at the master file. if (proc_rank == 0) then - open (unit=29, file=file_loc, form='formatted', position='rewind') - !write(29,*) 'lg_bubID, x, y, z, xPrev, yPrev, zPrev, xVel, yVel, ', & - ! 'zVel, radius, interfaceVelocity, equilibriumRadius', & - ! 'Rmax, Rmin, dphidt, pressure, mv, mg, betaT, betaC, time' - do i = 1, tot_data - id = int(MPI_IO_DATA_lg_bubbles(i, 1)) - inputvals(1:20) = MPI_IO_DATA_lg_bubbles(i, 2:21) - if (id > 0) then - write (29, 6) int(id), inputvals(1), inputvals(2), & - inputvals(3), inputvals(4), inputvals(5), inputvals(6), inputvals(7), & - inputvals(8), inputvals(9), inputvals(10), inputvals(11), & - inputvals(12), inputvals(13), inputvals(14), inputvals(15), & - inputvals(16), inputvals(17), inputvals(18), inputvals(19), & - inputvals(20), time_real -6 format(I6, 21(1x, E15.7)) - end if + + do i = 1, num_procs + write (meshnames(i), '(A,I0,A,I0,A)') '../p', i - 1, & + '/', t_step, '.silo:lag_bubbles' + meshtypes(i) = DB_POINTMESH end do - close (29) + err = DBSET2DSTRLEN(len(meshnames(1))) + err = DBPUTMMESH(dbroot, 'lag_bubbles', 16, & + num_procs, meshnames, & + len_trim(meshnames), & + meshtypes, DB_F77NULL, ierr) end if + err = DBPUTPM(dbfile, 'lag_bubbles', 11, 3, & + px, py, pz, nBub, & + DB_DOUBLE, DB_F77NULL, ierr) + + if (lag_id_wrt) call s_write_lag_variable_to_formatted_database_file('part_id', t_step, bub_id, nBub) + if (lag_vel_wrt) then + call s_write_lag_variable_to_formatted_database_file('part_vel1', t_step, vx, nBub) + call s_write_lag_variable_to_formatted_database_file('part_vel2', t_step, vy, nBub) + if (p > 0) call s_write_lag_variable_to_formatted_database_file('part_vel3', t_step, vz, nBub) + end if + if (lag_rad_wrt) call s_write_lag_variable_to_formatted_database_file('part_radius', t_step, radius, nBub) + if (lag_rvel_wrt) call s_write_lag_variable_to_formatted_database_file('part_rdot', t_step, rvel, nBub) + if (lag_r0_wrt) call s_write_lag_variable_to_formatted_database_file('part_r0', t_step, rnot, nBub) + if (lag_rmax_wrt) call s_write_lag_variable_to_formatted_database_file('part_rmax', t_step, rmax, nBub) + if (lag_rmin_wrt) call s_write_lag_variable_to_formatted_database_file('part_rmin', t_step, rmin, nBub) + if (lag_dphidt_wrt) call s_write_lag_variable_to_formatted_database_file('part_dphidt', t_step, dphidt, nBub) + if (lag_pres_wrt) call s_write_lag_variable_to_formatted_database_file('part_pressure', t_step, pressure, nBub) + if (lag_mv_wrt) call s_write_lag_variable_to_formatted_database_file('part_mv', t_step, mv, nBub) + if (lag_mg_wrt) call s_write_lag_variable_to_formatted_database_file('part_mg', t_step, mg, nBub) + if (lag_betaT_wrt) call s_write_lag_variable_to_formatted_database_file('part_betaT', t_step, betaT, nBub) + if (lag_betaC_wrt) call s_write_lag_variable_to_formatted_database_file('part_betaC', t_step, betaC, nBub) + + deallocate (bub_id, px, py, pz, ppx, ppy, ppz, vx, vy, vz, radius, & + rvel, rnot, rmax, rmin, dphidt, pressure, mv, mg, & + betaT, betaC) deallocate (MPI_IO_DATA_lg_bubbles) + else + call MPI_TYPE_CONTIGUOUS(0, mpi_p, view, ierr) + call MPI_TYPE_COMMIT(view, ierr) - end if + call MPI_FILE_OPEN(MPI_COMM_WORLD, file_loc, MPI_MODE_RDONLY, & + mpi_info_int, ifile, ierr) - call s_mpi_barrier() + ! Skip extended header + disp = int(sizeof(file_tot_part) + 2*sizeof(file_time) + sizeof(file_num_procs) + & + file_num_procs*sizeof(proc_bubble_counts(1)), MPI_OFFSET_KIND) + call MPI_FILE_SET_VIEW(ifile, disp, mpi_p, view, 'native', mpi_info_int, ierr) - call MPI_FILE_CLOSE(ifile, ierr) + call MPI_FILE_READ_ALL(ifile, dummy, 0, mpi_p, status, ierr) + + call MPI_FILE_CLOSE(ifile, ierr) + call MPI_TYPE_FREE(view, ierr) + + if (proc_rank == 0) then + + do i = 1, num_procs + write (meshnames(i), '(A,I0,A,I0,A)') '../p', i - 1, & + '/', t_step, '.silo:lag_bubbles' + meshtypes(i) = DB_POINTMESH + end do + err = DBSET2DSTRLEN(len(meshnames(1))) + err = DBPUTMMESH(dbroot, 'lag_bubbles', 16, & + num_procs, meshnames, & + len_trim(meshnames), & + meshtypes, DB_F77NULL, ierr) + end if + + err = DBSETEMPTYOK(1) + err = DBPUTPM(dbfile, 'lag_bubbles', 11, 3, & + dummy_data, dummy_data, dummy_data, 0, & + DB_DOUBLE, DB_F77NULL, ierr) + + if (lag_id_wrt) call s_write_lag_variable_to_formatted_database_file('part_id', t_step) + if (lag_vel_wrt) then + call s_write_lag_variable_to_formatted_database_file('part_vel1', t_step) + call s_write_lag_variable_to_formatted_database_file('part_vel2', t_step) + if (p > 0) call s_write_lag_variable_to_formatted_database_file('part_vel3', t_step) + end if + if (lag_rad_wrt) call s_write_lag_variable_to_formatted_database_file('part_radius', t_step) + if (lag_rvel_wrt) call s_write_lag_variable_to_formatted_database_file('part_rdot', t_step) + if (lag_r0_wrt) call s_write_lag_variable_to_formatted_database_file('part_r0', t_step) + if (lag_rmax_wrt) call s_write_lag_variable_to_formatted_database_file('part_rmax', t_step) + if (lag_rmin_wrt) call s_write_lag_variable_to_formatted_database_file('part_rmin', t_step) + if (lag_dphidt_wrt) call s_write_lag_variable_to_formatted_database_file('part_dphidt', t_step) + if (lag_pres_wrt) call s_write_lag_variable_to_formatted_database_file('part_pressure', t_step) + if (lag_mv_wrt) call s_write_lag_variable_to_formatted_database_file('part_mv', t_step) + if (lag_mg_wrt) call s_write_lag_variable_to_formatted_database_file('part_mg', t_step) + if (lag_betaT_wrt) call s_write_lag_variable_to_formatted_database_file('part_betaT', t_step) + if (lag_betaC_wrt) call s_write_lag_variable_to_formatted_database_file('part_betaC', t_step) + end if #endif - end subroutine s_write_lag_bubbles_results + end subroutine s_write_lag_bubbles_to_formatted_database_file + + subroutine s_write_lag_variable_to_formatted_database_file(varname, t_step, data, nBubs) + + character(len=*), intent(in) :: varname + integer, intent(in) :: t_step + real(wp), dimension(1:), intent(in), optional :: data + integer, intent(in), optional :: nBubs + + character(len=64), dimension(num_procs) :: var_names + integer, dimension(num_procs) :: var_types + real(wp) :: dummy_data + + integer :: ierr !< Generic flag used to identify and report database errors + integer :: i + + dummy_data = 0._wp + + if (present(nBubs) .and. present(data)) then + if (proc_rank == 0) then + do i = 1, num_procs + write (var_names(i), '(A,I0,A,I0,A)') '../p', i - 1, & + '/', t_step, '.silo:'//trim(varname) + var_types(i) = DB_POINTVAR + end do + err = DBSET2DSTRLEN(len(var_names(1))) + err = DBPUTMVAR(dbroot, trim(varname), len_trim(varname), & + num_procs, var_names, & + len_trim(var_names), & + var_types, DB_F77NULL, ierr) + end if + + err = DBPUTPV1(dbfile, trim(varname), len_trim(varname), & + 'lag_bubbles', 11, data, nBubs, DB_DOUBLE, DB_F77NULL, ierr) + else + if (proc_rank == 0) then + do i = 1, num_procs + write (var_names(i), '(A,I0,A,I0,A)') '../p', i - 1, & + '/', t_step, '.silo:'//trim(varname) + var_types(i) = DB_POINTVAR + end do + err = DBSET2DSTRLEN(len(var_names(1))) + err = DBSETEMPTYOK(1) + err = DBPUTMVAR(dbroot, trim(varname), len_trim(varname), & + num_procs, var_names, & + len_trim(var_names), & + var_types, DB_F77NULL, ierr) + end if + + err = DBSETEMPTYOK(1) + err = DBPUTPV1(dbfile, trim(varname), len_trim(varname), & + 'lag_bubbles', 11, dummy_data, 0, DB_DOUBLE, DB_F77NULL, ierr) + end if + + end subroutine s_write_lag_variable_to_formatted_database_file + impure subroutine s_write_intf_data_file(q_prim_vf) type(scalar_field), dimension(sys_size), intent(IN) :: q_prim_vf @@ -1443,7 +1734,7 @@ contains ! the offsets and the one bookkeeping the number of cell-boundaries ! in each active coordinate direction. Note that all these variables ! were only needed by Silo-HDF5 format for multidimensional data. - if (format == 1 .and. n > 0) then + if (format == 1) then deallocate (spatial_extents) deallocate (data_extents) deallocate (lo_offset) diff --git a/src/post_process/m_global_parameters.fpp b/src/post_process/m_global_parameters.fpp index 134d900f91..11c5415946 100644 --- a/src/post_process/m_global_parameters.fpp +++ b/src/post_process/m_global_parameters.fpp @@ -60,7 +60,6 @@ module m_global_parameters !> @name Cell-boundary locations in the x-, y- and z-coordinate directions !> @{ real(wp), allocatable, dimension(:) :: x_cb, x_root_cb, y_cb, z_cb - real(wp), allocatable, dimension(:) :: x_cb_s, y_cb_s, z_cb_s !> @} !> @name Cell-center locations in the x-, y- and z-coordinate directions @@ -178,6 +177,11 @@ module m_global_parameters integer, allocatable, dimension(:) :: proc_coords !< !! Processor coordinates in MPI_CART_COMM + type(int_bounds_info), dimension(3) :: nidx + + integer, allocatable, dimension(:, :, :) :: neighbor_ranks + !! Neighbor ranks for lagrangian particle communication + integer, allocatable, dimension(:) :: start_idx !< !! Starting cell-center index of local processor in global grid @@ -258,6 +262,24 @@ module m_global_parameters logical :: ib logical :: chem_wrt_Y(1:num_species) logical :: chem_wrt_T + logical :: lag_header + logical :: lag_txt_wrt + logical :: lag_db_wrt + logical :: lag_id_wrt + logical :: lag_pos_wrt + logical :: lag_pos_prev_wrt + logical :: lag_vel_wrt + logical :: lag_rad_wrt + logical :: lag_rvel_wrt + logical :: lag_r0_wrt + logical :: lag_rmax_wrt + logical :: lag_rmin_wrt + logical :: lag_dphidt_wrt + logical :: lag_pres_wrt + logical :: lag_mv_wrt + logical :: lag_mg_wrt + logical :: lag_betaT_wrt + logical :: lag_betaC_wrt !> @} real(wp), dimension(num_fluids_max) :: schlieren_alpha !< @@ -329,6 +351,8 @@ module m_global_parameters real(wp) :: Bx0 !< Constant magnetic field in the x-direction (1D) + real(wp) :: wall_time, wall_time_avg !< Wall time measurements + contains !> Assigns default values to user inputs prior to reading @@ -437,7 +461,43 @@ contains schlieren_wrt = .false. sim_data = .false. cf_wrt = .false. + lag_txt_wrt = .false. + lag_header = .true. + lag_db_wrt = .false. + lag_id_wrt = .true. + lag_pos_wrt = .true. + lag_pos_prev_wrt = .false. + lag_vel_wrt = .true. + lag_rad_wrt = .true. + lag_rvel_wrt = .false. + lag_r0_wrt = .false. + lag_rmax_wrt = .false. + lag_rmin_wrt = .false. + lag_dphidt_wrt = .false. + lag_pres_wrt = .false. + lag_mv_wrt = .false. + lag_mg_wrt = .false. + lag_betaT_wrt = .false. + lag_betaC_wrt = .false. ib = .false. + lag_txt_wrt = .false. + lag_header = .true. + lag_db_wrt = .false. + lag_id_wrt = .true. + lag_pos_wrt = .true. + lag_pos_prev_wrt = .false. + lag_vel_wrt = .true. + lag_rad_wrt = .true. + lag_rvel_wrt = .false. + lag_r0_wrt = .false. + lag_rmax_wrt = .false. + lag_rmin_wrt = .false. + lag_dphidt_wrt = .false. + lag_pres_wrt = .false. + lag_mv_wrt = .false. + lag_mg_wrt = .false. + lag_betaT_wrt = .false. + lag_betaC_wrt = .false. schlieren_alpha = dflt_real @@ -812,7 +872,7 @@ contains ! in the Silo-HDF5 format. If this is the case, one must also verify ! whether the raw simulation data is 2D or 3D. In the 2D case, size ! of the z-coordinate direction ghost zone layer must be zeroed out. - if (num_procs == 1 .or. format /= 1 .or. n == 0) then + if (num_procs == 1 .or. format /= 1) then offset_x%beg = 0 offset_x%end = 0 @@ -821,6 +881,13 @@ contains offset_z%beg = 0 offset_z%end = 0 + elseif (n == 0) then + + offset_y%beg = 0 + offset_y%end = 0 + offset_z%beg = 0 + offset_z%end = 0 + elseif (p == 0) then offset_z%beg = 0 @@ -855,17 +922,7 @@ contains idwbuff(3)%end = idwint(3)%end - idwbuff(3)%beg ! Allocating single precision grid variables if needed - if (precision == 1) then - allocate (x_cb_s(-1 - offset_x%beg:m + offset_x%end)) - if (n > 0) then - allocate (y_cb_s(-1 - offset_y%beg:n + offset_y%end)) - if (p > 0) then - allocate (z_cb_s(-1 - offset_z%beg:p + offset_z%end)) - end if - end if - else - allocate (x_cc_s(-buff_size:m + buff_size)) - end if + allocate (x_cc_s(-buff_size:m + buff_size)) ! Allocating the grid variables in the x-coordinate direction allocate (x_cb(-1 - offset_x%beg:m + offset_x%end)) diff --git a/src/post_process/m_mpi_proxy.fpp b/src/post_process/m_mpi_proxy.fpp index 8699b97fe4..795936250e 100644 --- a/src/post_process/m_mpi_proxy.fpp +++ b/src/post_process/m_mpi_proxy.fpp @@ -106,10 +106,19 @@ contains & 'surface_tension', 'hyperelasticity', 'bubbles_lagrange', & & 'output_partial_domain', 'relativity', 'cont_damage', 'bc_io', & & 'down_sample' ] - call MPI_BCAST(${VAR}$, 1, MPI_LOGICAL, 0, MPI_COMM_WORLD, ierr) #:endfor + if (bubbles_lagrange) then + #:for VAR in ['lag_header', 'lag_txt_wrt', 'lag_db_wrt', 'lag_id_wrt', & + & 'lag_pos_wrt', 'lag_pos_prev_wrt', 'lag_vel_wrt', 'lag_rad_wrt', & + & 'lag_rvel_wrt', 'lag_r0_wrt', 'lag_rmax_wrt', 'lag_rmin_wrt', & + & 'lag_dphidt_wrt', 'lag_pres_wrt', 'lag_mv_wrt', 'lag_mg_wrt', & + & 'lag_betaT_wrt', 'lag_betaC_wrt', 'bc_io', 'down_sample' ] + call MPI_BCAST(${VAR}$, 1, MPI_LOGICAL, 0, MPI_COMM_WORLD, ierr) + #:endfor + end if + call MPI_BCAST(flux_wrt(1), 3, MPI_LOGICAL, 0, MPI_COMM_WORLD, ierr) call MPI_BCAST(omega_wrt(1), 3, MPI_LOGICAL, 0, MPI_COMM_WORLD, ierr) call MPI_BCAST(mom_wrt(1), 3, MPI_LOGICAL, 0, MPI_COMM_WORLD, ierr) @@ -226,7 +235,7 @@ contains ierr) end if ! Simulation is 2D - else + elseif (n > 0) then ! Minimum spatial extent in the x-direction call MPI_GATHERV(minval(x_cb), 1, mpi_p, & @@ -251,7 +260,20 @@ contains spatial_extents(4, 0), recvcounts, 4*displs, & mpi_p, 0, MPI_COMM_WORLD, & ierr) + ! Simulation is 1D + else + + ! Minimum spatial extent in the x-direction + call MPI_GATHERV(minval(x_cb), 1, mpi_p, & + spatial_extents(1, 0), recvcounts, 4*displs, & + mpi_p, 0, MPI_COMM_WORLD, & + ierr) + ! Maximum spatial extent in the x-direction + call MPI_GATHERV(maxval(x_cb), 1, mpi_p, & + spatial_extents(2, 0), recvcounts, 4*displs, & + mpi_p, 0, MPI_COMM_WORLD, & + ierr) end if #endif diff --git a/src/post_process/m_start_up.f90 b/src/post_process/m_start_up.f90 index c5e7fe31be..89e4d838eb 100644 --- a/src/post_process/m_start_up.f90 +++ b/src/post_process/m_start_up.f90 @@ -90,7 +90,11 @@ impure subroutine s_read_input_file cfl_target, surface_tension, bubbles_lagrange, & sim_data, hyperelasticity, Bx0, relativity, cont_damage, & num_bc_patches, igr, igr_order, down_sample, recon_type, & - muscl_order + muscl_order, lag_header, lag_txt_wrt, lag_db_wrt, & + lag_id_wrt, lag_pos_wrt, lag_pos_prev_wrt, lag_vel_wrt, & + lag_rad_wrt, lag_rvel_wrt, lag_r0_wrt, lag_rmax_wrt, & + lag_rmin_wrt, lag_dphidt_wrt, lag_pres_wrt, lag_mv_wrt, & + lag_mg_wrt, lag_betaT_wrt, lag_betaC_wrt ! Inquiring the status of the post_process.inp file file_loc = 'post_process.inp' @@ -176,15 +180,15 @@ impure subroutine s_perform_time_step(t_step) integer, intent(inout) :: t_step if (proc_rank == 0) then if (cfl_dt) then - print '(" [", I3, "%] Saving ", I8, " of ", I0, "")', & + print '(" [", I3, "%] Saving ", I8, " of ", I0, " Time Avg = ", ES16.6, " Time/step = ", ES12.6, "")', & int(ceiling(100._wp*(real(t_step - n_start)/(n_save)))), & - t_step, n_save + t_step, n_save, wall_time_avg, wall_time else - print '(" [", I3, "%] Saving ", I8, " of ", I0, " @ t_step = ", I0, "")', & + print '(" [", I3, "%] Saving ", I8, " of ", I0, " @ t_step = ", I8, " Time Avg = ", ES16.6, " Time/step = ", ES12.6, "")', & int(ceiling(100._wp*(real(t_step - t_step_start)/(t_step_stop - t_step_start + 1)))), & (t_step - t_step_start)/t_step_save + 1, & (t_step_stop - t_step_start)/t_step_save + 1, & - t_step + t_step, wall_time_avg, wall_time end if end if @@ -719,7 +723,8 @@ impure subroutine s_save_data(t_step, varname, pres, c, H) call s_write_variable_to_formatted_database_file(varname, t_step) varname(:) = ' ' - call s_write_lag_bubbles_results(t_step) !! Individual bubble evolution + if (lag_txt_wrt) call s_write_lag_bubbles_results_to_text(t_step) ! text output + if (lag_db_wrt) call s_write_lag_bubbles_to_formatted_database_file(t_step) ! silo file output end if if (sim_data .and. proc_rank == 0) then @@ -736,9 +741,7 @@ impure subroutine s_initialize_modules ! Computation of parameters, allocation procedures, and/or any other tasks ! needed to properly setup the modules call s_initialize_global_parameters_module() - if (bubbles_euler .and. nb > 1) then - call s_simpson(weight, R0) - end if + if (bubbles_euler .and. nb > 1) call s_simpson(weight, R0) if (bubbles_euler .and. .not. polytropic) then call s_initialize_nonpoly() end if diff --git a/src/post_process/p_main.fpp b/src/post_process/p_main.fpp index 59b69e3b07..12d6353c67 100644 --- a/src/post_process/p_main.fpp +++ b/src/post_process/p_main.fpp @@ -26,6 +26,7 @@ program p_main real(wp) :: pres real(wp) :: c real(wp) :: H + real(wp) :: start, finish call s_initialize_mpi_domain() @@ -49,10 +50,22 @@ program p_main ! available step. To avoid this, we force synchronization here. call s_mpi_barrier() + call cpu_time(start) + call s_perform_time_step(t_step) call s_save_data(t_step, varname, pres, c, H) + call cpu_time(finish) + + wall_time = abs(finish - start) + + if (t_step >= 2) then + wall_time_avg = (wall_time + (t_step - 2)*wall_time_avg)/(t_step - 1) + else + wall_time_avg = 0._wp + end if + if (cfl_dt) then if (t_step == n_save - 1) then exit diff --git a/src/pre_process/m_assign_variables.fpp b/src/pre_process/m_assign_variables.fpp index f60573b1b9..ca0c7c1c16 100644 --- a/src/pre_process/m_assign_variables.fpp +++ b/src/pre_process/m_assign_variables.fpp @@ -684,7 +684,7 @@ contains if (surface_tension) then q_prim_vf(c_idx)%sf(j, k, l) = eta*patch_icpp(patch_id)%cf_val + & - (1._wp - eta)*patch_icpp(smooth_patch_id)%cf_val + (1._wp - eta)*orig_prim_vf(c_idx) end if ! Updating the patch identities bookkeeping variable diff --git a/src/pre_process/m_global_parameters.fpp b/src/pre_process/m_global_parameters.fpp index 79c23f2917..e2d0a2b41c 100644 --- a/src/pre_process/m_global_parameters.fpp +++ b/src/pre_process/m_global_parameters.fpp @@ -28,6 +28,9 @@ module m_global_parameters logical :: old_ic, non_axis_sym !< Use existing IC data integer :: t_step_old, t_step_start !< Existing IC/grid folder + character(LEN=path_len) :: interface_file + real(wp) :: normFac, normMag, g0, p0 + logical :: cfl_adap_dt, cfl_const_dt, cfl_dt integer :: n_start, n_start_old @@ -131,6 +134,21 @@ module m_global_parameters ! Stands for "InDices With BUFFer". type(int_bounds_info) :: idwbuff(1:3) + integer :: fd_order !< + !! The order of the finite-difference (fd) approximations of the first-order + !! derivatives that need to be evaluated when the CoM or flow probe data + !! files are to be written at each time step + + integer :: fd_number !< + !! The finite-difference number is given by MAX(1, fd_order/2). Essentially, + !! it is a measure of the half-size of the finite-difference stencil for the + !! selected order of accuracy. + + !> @name lagrangian subgrid bubble parameters + !> @{! + type(bubbles_lagrange_parameters) :: lag_params !< Lagrange bubbles' parameters + !> @} + type(int_bounds_info) :: bc_x, bc_y, bc_z !< !! Boundary conditions in the x-, y- and z-coordinate directions @@ -175,6 +193,11 @@ module m_global_parameters integer, allocatable, dimension(:) :: proc_coords !< !! Processor coordinates in MPI_CART_COMM + type(int_bounds_info), dimension(3) :: nidx + + integer, allocatable, dimension(:, :, :) :: neighbor_ranks + !! Neighbor ranks for lagrangian particle communication + integer, allocatable, dimension(:) :: start_idx !< !! Starting cell-center index of local processor in global grid @@ -303,6 +326,11 @@ contains old_ic = .false. t_step_old = dflt_int t_step_start = dflt_int + interface_file = '.' + normFac = dflt_real + normMag = dflt_real + g0 = dflt_real + p0 = dflt_real cfl_adap_dt = .false. cfl_const_dt = .false. @@ -398,6 +426,28 @@ contains ! Initial condition parameters num_patches = dflt_int + fd_order = dflt_int + lag_params%cluster_type = dflt_int + lag_params%pressure_corrector = .false. + lag_params%smooth_type = dflt_int + lag_params%heatTransfer_model = .false. + lag_params%massTransfer_model = .false. + lag_params%write_bubbles = .false. + lag_params%write_bubbles_stats = .false. + lag_params%nBubs_glb = dflt_int + lag_params%vel_model = dflt_int + lag_params%drag_model = dflt_int + lag_params%c_d = dflt_real + lag_params%epsilonb = 1._wp + lag_params%charwidth = dflt_real + lag_params%valmaxvoid = dflt_real + lag_params%c0 = dflt_real + lag_params%rho0 = dflt_real + lag_params%T0 = dflt_real + lag_params%Thost = dflt_real + lag_params%x0 = dflt_real + lag_params%diffcoefvap = dflt_real + do i = 1, num_patches_max patch_icpp(i)%geometry = dflt_int patch_icpp(i)%model_scale(:) = 1._wp @@ -893,11 +943,15 @@ contains chemxb = species_idx%beg chemxe = species_idx%end + if (lag_params%vel_model /= 0) then + fd_number = max(1, fd_order/2) + end if + call s_configure_coordinate_bounds(recon_type, weno_polyn, muscl_polyn, & igr_order, buff_size, & idwint, idwbuff, viscous, & bubbles_lagrange, m, n, p, & - num_dims, igr, ib) + num_dims, igr, ib, fd_number) #ifdef MFC_MPI diff --git a/src/pre_process/m_mpi_proxy.fpp b/src/pre_process/m_mpi_proxy.fpp index c4ab3aa58b..21c05f6eed 100644 --- a/src/pre_process/m_mpi_proxy.fpp +++ b/src/pre_process/m_mpi_proxy.fpp @@ -40,14 +40,15 @@ contains ! Logistics call MPI_BCAST(case_dir, len(case_dir), MPI_CHARACTER, 0, MPI_COMM_WORLD, ierr) + call MPI_BCAST(interface_file, len(interface_file), MPI_CHARACTER, 0, MPI_COMM_WORLD, ierr) #:for VAR in ['t_step_old', 't_step_start', 'm', 'n', 'p', 'm_glb', 'n_glb', 'p_glb', & & 'loops_x', 'loops_y', 'loops_z', 'model_eqns', 'num_fluids', & & 'weno_order', 'precision', 'perturb_flow_fluid', & & 'perturb_sph_fluid', 'num_patches', 'thermal', 'nb', 'dist_type',& & 'relax_model', 'num_ibs', 'n_start', 'elliptic_smoothing_iters', & - & 'num_bc_patches', 'mixlayer_perturb_nk', 'recon_type', & - & 'muscl_order', 'igr_order' ] + & 'num_bc_patches', 'mixlayer_perturb_nk', 'fd_order', & + & 'lag_params%vel_model', 'recon_type', 'muscl_order', 'igr_order' ] call MPI_BCAST(${VAR}$, 1, MPI_INTEGER, 0, MPI_COMM_WORLD, ierr) #:endfor @@ -71,7 +72,7 @@ contains & 'perturb_flow_mag', 'pref', 'rhoref', 'poly_sigma', 'R0ref', & & 'Web', 'Ca', 'Re_inv', 'sigR', 'sigV', 'rhoRV', 'palpha_eps', & & 'ptgalpha_eps', 'sigma', 'pi_fac', 'mixlayer_vel_coef', 'Bx0', & - & 'mixlayer_perturb_k0'] + & 'mixlayer_perturb_k0', 'normFac', 'normMag', 'p0', 'g0'] call MPI_BCAST(${VAR}$, 1, mpi_p, 0, MPI_COMM_WORLD, ierr) #:endfor diff --git a/src/pre_process/m_start_up.fpp b/src/pre_process/m_start_up.fpp index 85b3d18bd5..19c004b7cf 100644 --- a/src/pre_process/m_start_up.fpp +++ b/src/pre_process/m_start_up.fpp @@ -151,8 +151,9 @@ contains n_start_old, surface_tension, hyperelasticity, pre_stress, & elliptic_smoothing, elliptic_smoothing_iters, & viscous, bubbles_lagrange, bc_x, bc_y, bc_z, num_bc_patches, & - patch_bc, Bx0, relativity, cont_damage, igr, igr_order, & - down_sample, recon_type, muscl_order + patch_bc, Bx0, relativity, cont_damage, fd_order, lag_params, & + interface_file, normFac, normMag, g0, p0, igr, igr_order, & + igr, igr_order, down_sample, recon_type, muscl_order ! Inquiring the status of the pre_process.inp file file_loc = 'pre_process.inp' diff --git a/src/simulation/m_body_forces.fpp b/src/simulation/m_body_forces.fpp index 41a103fcee..26e2663164 100644 --- a/src/simulation/m_body_forces.fpp +++ b/src/simulation/m_body_forces.fpp @@ -56,16 +56,13 @@ contains subroutine s_compute_acceleration(t) real(wp), intent(in) :: t + accel_bf(:) = 0._wp - if (m > 0) then - accel_bf(1) = g_x + k_x*sin(w_x*t - p_x) - if (n > 0) then - accel_bf(2) = g_y + k_y*sin(w_y*t - p_y) - if (p > 0) then - accel_bf(3) = g_z + k_z*sin(w_z*t - p_z) - end if + #:for DIR, XYZ in [(1, 'x'), (2, 'y'), (3, 'z')] + if (bf_${XYZ}$) then + accel_bf(${DIR}$) = g_${XYZ}$+k_${XYZ}$*sin(w_${XYZ}$*t - p_${XYZ}$) end if - end if + #:endfor $:GPU_UPDATE(device='[accel_bf]') diff --git a/src/simulation/m_bubbles.fpp b/src/simulation/m_bubbles.fpp index 0ec758dc22..2553f4d29c 100644 --- a/src/simulation/m_bubbles.fpp +++ b/src/simulation/m_bubbles.fpp @@ -17,6 +17,8 @@ module m_bubbles use m_helper_basic !< Functions to compare floating point numbers + use m_bubbles_EL_kernels + implicit none real(wp) :: chi_vw !< Bubble wall properties (Ando 2010) @@ -49,7 +51,10 @@ contains real(wp) :: fCpbw, fCpinf, fCpinf_dot, fH, fHdot, c_gas, c_liquid real(wp) :: f_rddot - if (bubble_model == 1) then + if (bubble_model == 0) then + ! Particle + f_rddot = 0._wp + else if (bubble_model == 1) then ! Gilmore bubbles fCpinf = fP - pref fCpbw = f_cpbw(fR0, fR, fV, fpb) @@ -72,6 +77,9 @@ contains ! Rayleigh-Plesset bubbles fCpbw = f_cpbw_KM(fR0, fR, fV, fpb) f_rddot = f_rddot_RP(fP, fRho, fR, fV, fCpbw) + else + ! Default: No bubble dynamics + f_rddot = 0._wp end if end function f_rddot @@ -445,7 +453,7 @@ contains !! @param fRho Current density !! @param fP Current driving pressure !! @param fR Current bubble radius - !! @param fV Current bubble velocity + !! @param fV Current bubble radial velocity !! @param fR0 Equilibrium bubble radius !! @param fpb Internal bubble pressure !! @param fpbdot Time-derivative of internal bubble pressure @@ -464,7 +472,8 @@ contains pure subroutine s_advance_step(fRho, fP, fR, fV, fR0, fpb, fpbdot, alf, & fntait, fBtait, f_bub_adv_src, f_divu, & bub_id, fmass_v, fmass_n, fbeta_c, & - fbeta_t, fCson, adap_dt_stop) + fbeta_t, fCson, adap_dt_stop, fRe, fPos, & + fVel, cell, q_prim_vf) $:GPU_ROUTINE(function_name='s_advance_step',parallelism='[seq]', & & cray_inline=True) @@ -474,14 +483,18 @@ contains integer, intent(in) :: bub_id real(wp), intent(in) :: fmass_n, fbeta_c, fbeta_t, fCson integer, intent(inout) :: adap_dt_stop + real(wp), intent(inout), dimension(3), optional :: fPos, fVel + real(wp), intent(in), optional :: fRe + integer, intent(in), dimension(3), optional :: cell + type(scalar_field), intent(in), dimension(sys_size), optional :: q_prim_vf real(wp), dimension(5) :: err !< Error estimates for adaptive time stepping real(wp) :: t_new !< Updated time step size real(wp) :: h !< Time step size real(wp), dimension(4) :: myR_tmp1, myV_tmp1, myR_tmp2, myV_tmp2 !< Bubble radius, radial velocity, and radial acceleration for the inner loop real(wp), dimension(4) :: myPb_tmp1, myMv_tmp1, myPb_tmp2, myMv_tmp2 !< Gas pressure and vapor mass for the inner loop (EL) - real(wp) :: fR2, fV2, fpb2, fmass_v2 - integer :: iter_count + real(wp) :: fR2, fV2, fpb2, fmass_v2, vTemp, aTemp, f_bTemp + integer :: l, iter_count call s_initial_substep_h(fRho, fP, fR, fV, fR0, fpb, fpbdot, alf, & fntait, fBtait, f_bub_adv_src, f_divu, fCson, h) @@ -552,6 +565,30 @@ contains ! Update pb and mass_v fpb = myPb_tmp1(4) fmass_v = myMv_tmp1(4) + + if (lag_params%vel_model == 1) then + do l = 1, num_dims + vTemp = f_interpolate_velocity(fR, cell, l, q_prim_vf) + fPos(l) = fPos(l) + h*vTemp + fVel(l) = vTemp + end do + elseif (lag_params%vel_model == 2) then + do l = 1, num_dims + f_bTemp = f_get_bubble_force(fPos(l), fR, fV, fVel, fmass_n, fmass_v, & + fRe, fRho, cell, l, q_prim_vf) + aTemp = f_bTemp/(fmass_n + fmass_v) + fPos(l) = fPos(l) + h * fVel(l) + fVel(l) = fVel(l) + h * aTemp + end do + elseif (lag_params%vel_model == 3) then + do l = 1, num_dims + f_bTemp = f_get_bubble_force(fPos(l), fR, fV, fVel, fmass_n, fmass_v, & + fRe, fRho, cell, l, q_prim_vf) + aTemp = 2._wp * f_bTemp / (fmass_n + fmass_v) - 3 * fV * fVel(l) / fR + fPos(l) = fPos(l) + h * fVel(l) + fVel(l) = fVel(l) + h * aTemp + end do + end if end if ! Update step size for the next sub-step diff --git a/src/simulation/m_bubbles_EL.fpp b/src/simulation/m_bubbles_EL.fpp index f5dbc1c628..35d70eb9f6 100644 --- a/src/simulation/m_bubbles_EL.fpp +++ b/src/simulation/m_bubbles_EL.fpp @@ -11,6 +11,8 @@ module m_bubbles_EL use m_mpi_proxy !< Message passing interface (MPI) module proxy + use m_mpi_common + use m_bubbles_EL_kernels !< Definitions of the kernel functions use m_bubbles !< General bubble dynamics procedures @@ -27,6 +29,8 @@ module m_bubbles_EL use m_helper + use m_ibm + implicit none !(nBub) @@ -74,6 +78,12 @@ module m_bubbles_EL $:GPU_DECLARE(create='[nBubs,Rmax_glb,Rmin_glb,q_beta,q_beta_idx]') + integer, allocatable, dimension(:) :: keep_bubble, prefix_sum + integer, allocatable, dimension(:,:) :: wrap_bubble_loc, wrap_bubble_dir + integer :: active_bubs + $:GPU_DECLARE(create='[keep_bubble, prefix_sum, active_bubs]') + $:GPU_DECLARE(create='[wrap_bubble_loc, wrap_bubble_dir]') + contains !> Initializes the lagrangian subgrid bubble solver @@ -102,6 +112,20 @@ contains call s_mpi_abort('Please check the lag_params%solver_approach input') end if + pcomm_coords(1)%beg = x_cb(buff_size - fd_number - 1) + pcomm_coords(1)%end = x_cb(m - buff_size + fd_number) + $:GPU_UPDATE(device='[pcomm_coords(1)]') + if (n > 0) then + pcomm_coords(2)%beg = y_cb(buff_size - fd_number - 1) + pcomm_coords(2)%end = y_cb(n - buff_size + fd_number) + $:GPU_UPDATE(device='[pcomm_coords(2)]') + if (p > 0) then + pcomm_coords(3)%beg = z_cb(buff_size - fd_number - 1) + pcomm_coords(3)%end = z_cb(p - buff_size + fd_number) + $:GPU_UPDATE(device='[pcomm_coords(3)]') + end if + end if + $:GPU_UPDATE(device='[lag_num_ts, q_beta_idx]') @:ALLOCATE(q_beta%vf(1:q_beta_idx)) @@ -117,7 +141,6 @@ contains ! Allocating space for lagrangian variables nBubs_glb = lag_params%nBubs_glb - @:ALLOCATE(lag_id(1:nBubs_glb, 1:2)) @:ALLOCATE(bub_R0(1:nBubs_glb)) @:ALLOCATE(Rmax_stats(1:nBubs_glb)) @:ALLOCATE(Rmin_stats(1:nBubs_glb)) @@ -125,6 +148,7 @@ contains @:ALLOCATE(gas_betaT(1:nBubs_glb)) @:ALLOCATE(gas_betaC(1:nBubs_glb)) @:ALLOCATE(bub_dphidt(1:nBubs_glb)) + @:ALLOCATE(lag_id(1:nBubs_glb, 1:2)) @:ALLOCATE(gas_p(1:nBubs_glb, 1:2)) @:ALLOCATE(gas_mv(1:nBubs_glb, 1:2)) @:ALLOCATE(intfc_rad(1:nBubs_glb, 1:2)) @@ -140,8 +164,13 @@ contains @:ALLOCATE(mtn_dposdt(1:nBubs_glb, 1:3, 1:lag_num_ts)) @:ALLOCATE(mtn_dveldt(1:nBubs_glb, 1:3, 1:lag_num_ts)) + @:ALLOCATE(keep_bubble(1:nBubs_glb), prefix_sum(1:nBubs_glb)) + @:ALLOCATE(wrap_bubble_loc(1:nBubs_glb, 1:num_dims), wrap_bubble_dir(1:nBubs_glb, 1:num_dims)) + if (adap_dt .and. f_is_default(adap_dt_tol)) adap_dt_tol = dflt_adap_dt_tol + if (num_procs > 1) call s_initialize_particles_mpi(lag_num_ts) + ! Starting bubbles call s_start_lagrange_inputs() call s_read_input_bubbles(q_cons_vf) @@ -226,7 +255,7 @@ contains do while (ios == 0) read (94, *, iostat=ios) (inputBubble(i), i=1, 8) if (ios /= 0) cycle - indomain = particle_in_domain(inputBubble(1:3)) + indomain = particle_in_domain_physical(inputBubble(1:3)) id = id + 1 if (id > lag_params%nBubs_glb .and. proc_rank == 0) then call s_mpi_abort("Current number of bubbles is larger than nBubs_glb") @@ -250,6 +279,22 @@ contains print *, " Lagrange bubbles running, in proc", proc_rank, "number:", bub_id, "/", id + call s_mpi_reduce_int_sum(bub_id, bub_id) + + if (proc_rank == 0) then + if (bub_id == 0) call s_mpi_abort('No bubbles in the domain. Check input/lag_bubbles.dat') + end if + + if (num_procs > 1) then + call s_add_particles_to_transfer_list(nBubs, mtn_pos(:, :, 1)) + call s_mpi_sendrecv_particles(bub_R0, Rmax_stats, Rmin_stats, gas_mg, gas_betaT, & + gas_betaC, bub_dphidt, lag_id, gas_p, gas_mv, & + intfc_rad, intfc_vel, mtn_pos, mtn_posPrev, mtn_vel, & + mtn_s, intfc_draddt, intfc_dveldt, gas_dpdt, & + gas_dmvdt, mtn_dposdt, mtn_dveldt, lag_num_ts, nBubs, & + dest=1) + end if + $:GPU_UPDATE(device='[bubbles_lagrange, lag_params]') $:GPU_UPDATE(device='[lag_id,bub_R0,Rmax_stats,Rmin_stats,gas_mg, & @@ -321,28 +366,24 @@ contains mtn_posPrev(bub_id, 1:3, 1) = mtn_pos(bub_id, 1:3, 1) end if - cell = -buff_size + cell = fd_number - buff_size call s_locate_cell(mtn_pos(bub_id, 1:3, 1), cell, mtn_s(bub_id, 1:3, 1)) - ! Check if the bubble is located in the ghost cell of a symmetric boundary - if ((bc_x%beg == BC_REFLECTIVE .and. cell(1) < 0) .or. & - (bc_x%end == BC_REFLECTIVE .and. cell(1) > m) .or. & - (bc_y%beg == BC_REFLECTIVE .and. cell(2) < 0) .or. & - (bc_y%end == BC_REFLECTIVE .and. cell(2) > n)) then - call s_mpi_abort("Lagrange bubble is in the ghost cells of a symmetric boundary.") + ! Check if the bubble is located in the ghost cell of a symmetric, or wall boundary + if ((any(bc_x%beg == (/BC_REFLECTIVE, BC_CHAR_SLIP_WALL, BC_SLIP_WALL, BC_NO_SLIP_WALL/)) .and. cell(1) < 0) .or. & + (any(bc_x%end == (/BC_REFLECTIVE, BC_CHAR_SLIP_WALL, BC_SLIP_WALL, BC_NO_SLIP_WALL/)) .and. cell(1) > m) .or. & + (any(bc_y%beg == (/BC_REFLECTIVE, BC_CHAR_SLIP_WALL, BC_SLIP_WALL, BC_NO_SLIP_WALL/)) .and. cell(2) < 0) .or. & + (any(bc_y%end == (/BC_REFLECTIVE, BC_CHAR_SLIP_WALL, BC_SLIP_WALL, BC_NO_SLIP_WALL/)) .and. cell(2) > n)) then + call s_mpi_abort("Lagrange bubble is in the ghost cells of a symmetric or wall boundary.") end if if (p > 0) then - if ((bc_z%beg == BC_REFLECTIVE .and. cell(3) < 0) .or. & - (bc_z%end == BC_REFLECTIVE .and. cell(3) > p)) then - call s_mpi_abort("Lagrange bubble is in the ghost cells of a symmetric boundary.") + if ((any(bc_z%beg == (/BC_REFLECTIVE, BC_CHAR_SLIP_WALL, BC_SLIP_WALL, BC_NO_SLIP_WALL/)) .and. cell(3) < 0) .or. & + (any(bc_z%end == (/BC_REFLECTIVE, BC_CHAR_SLIP_WALL, BC_SLIP_WALL, BC_NO_SLIP_WALL/)) .and. cell(3) > p)) then + call s_mpi_abort("Lagrange bubble is in the ghost cells of a symmetric or wall boundary.") end if end if - ! If particle is in the ghost cells, find the closest non-ghost cell - cell(1) = min(max(cell(1), 0), m) - cell(2) = min(max(cell(2), 0), n) - if (p > 0) cell(3) = min(max(cell(3), 0), p) call s_convert_to_mixture_variables(q_cons_vf, cell(1), cell(2), cell(3), & rhol, gamma, pi_inf, qv, Re) dynP = 0._wp @@ -405,6 +446,8 @@ contains integer, intent(inout) :: bub_id, save_count character(LEN=path_len + 2*name_len) :: file_loc + real(wp) :: file_time, file_dt + integer :: file_num_procs, file_tot_part, tot_part #ifdef MFC_MPI real(wp), dimension(20) :: inputvals @@ -419,81 +462,146 @@ contains integer :: ifile, ierr, tot_data, id integer :: i - write (file_loc, '(a,i0,a)') 'lag_bubbles_mpi_io_', save_count, '.dat' + integer, dimension(:), allocatable :: proc_bubble_counts + real(wp), dimension(1:1, 1:lag_io_vars) :: dummy + dummy = 0._wp + + ! Construct file path + write (file_loc, '(A,I0,A)') 'lag_bubbles_', save_count, '.dat' file_loc = trim(case_dir)//'/restart_data'//trim(mpiiofs)//trim(file_loc) - inquire (file=trim(file_loc), exist=file_exist) - - if (file_exist) then - if (proc_rank == 0) then - open (9, file=trim(file_loc), form='unformatted', status='unknown') - read (9) tot_data, mytime, dt - close (9) - print *, 'Reading lag_bubbles_mpi_io: ', tot_data, mytime, dt - end if - else - print '(a)', trim(file_loc)//' is missing. exiting.' - call s_mpi_abort + + ! Check if file exists + inquire (FILE=trim(file_loc), EXIST=file_exist) + if (.not. file_exist) then + call s_mpi_abort('Restart file '//trim(file_loc)//' does not exist!') end if - call MPI_BCAST(tot_data, 1, MPI_INTEGER, 0, MPI_COMM_WORLD, ierr) - call MPI_BCAST(mytime, 1, mpi_p, 0, MPI_COMM_WORLD, ierr) - call MPI_BCAST(dt, 1, mpi_p, 0, MPI_COMM_WORLD, ierr) + if (.not. parallel_io) return + + if (proc_rank == 0) then + call MPI_FILE_OPEN(MPI_COMM_SELF, file_loc, MPI_MODE_RDONLY, & + mpi_info_int, ifile, ierr) + + call MPI_FILE_READ(ifile, file_tot_part, 1, MPI_INTEGER, status, ierr) + call MPI_FILE_READ(ifile, file_time, 1, mpi_p, status, ierr) + call MPI_FILE_READ(ifile, file_dt, 1, mpi_p, status, ierr) + call MPI_FILE_READ(ifile, file_num_procs, 1, MPI_INTEGER, status, ierr) + + call MPI_FILE_CLOSE(ifile, ierr) + end if + + call MPI_BCAST(file_tot_part, 1, MPI_INTEGER, 0, MPI_COMM_WORLD, ierr) + call MPI_BCAST(file_time, 1, mpi_p, 0, MPI_COMM_WORLD, ierr) + call MPI_BCAST(file_dt, 1, mpi_p, 0, MPI_COMM_WORLD, ierr) + call MPI_BCAST(file_num_procs, 1, MPI_INTEGER, 0, MPI_COMM_WORLD, ierr) + + allocate (proc_bubble_counts(file_num_procs)) + + if (proc_rank == 0) then + call MPI_FILE_OPEN(MPI_COMM_SELF, file_loc, MPI_MODE_RDONLY, & + mpi_info_int, ifile, ierr) + + ! Skip to processor counts position + disp = int(sizeof(file_tot_part) + 2*sizeof(file_time) + sizeof(file_num_procs), & + MPI_OFFSET_KIND) + call MPI_FILE_SEEK(ifile, disp, MPI_SEEK_SET, ierr) + call MPI_FILE_READ(ifile, proc_bubble_counts, file_num_procs, MPI_INTEGER, status, ierr) + + call MPI_FILE_CLOSE(ifile, ierr) + end if + + call MPI_BCAST(proc_bubble_counts, file_num_procs, MPI_INTEGER, 0, MPI_COMM_WORLD, ierr) + + ! Set time variables from file + mytime = file_time + dt = file_dt + + bub_id = proc_bubble_counts(proc_rank + 1) - gsizes(1) = tot_data - gsizes(2) = 21 - lsizes(1) = tot_data - lsizes(2) = 21 start_idx_part(1) = 0 + do i = 1, proc_rank + start_idx_part(1) = start_idx_part(1) + proc_bubble_counts(i) + end do + start_idx_part(2) = 0 + lsizes(1) = bub_id + lsizes(2) = lag_io_vars - call MPI_type_CREATE_SUBARRAY(2, gsizes, lsizes, start_idx_part, & - MPI_ORDER_FORTRAN, mpi_p, view, ierr) - call MPI_type_COMMIT(view, ierr) + gsizes(1) = file_tot_part + gsizes(2) = lag_io_vars - ! Open the file to write all flow variables - write (file_loc, '(a,i0,a)') 'lag_bubbles_', save_count, '.dat' - file_loc = trim(case_dir)//'/restart_data'//trim(mpiiofs)//trim(file_loc) - inquire (file=trim(file_loc), exist=particle_file) + if (bub_id > 0) then + + allocate (MPI_IO_DATA_lag_bubbles(bub_id, 1:lag_io_vars)) - if (particle_file) then - call MPI_FILE_open(MPI_COMM_WORLD, file_loc, MPI_MODE_RDONLY, & + call MPI_TYPE_CREATE_SUBARRAY(2, gsizes, lsizes, start_idx_part, & + MPI_ORDER_FORTRAN, mpi_p, view, ierr) + call MPI_TYPE_COMMIT(view, ierr) + + call MPI_FILE_OPEN(MPI_COMM_WORLD, file_loc, MPI_MODE_RDONLY, & mpi_info_int, ifile, ierr) - disp = 0._wp - call MPI_FILE_SET_VIEW(ifile, disp, mpi_p, view, & - 'native', mpi_info_null, ierr) - allocate (MPI_IO_DATA_lag_bubbles(tot_data, 1:21)) - call MPI_FILE_read_ALL(ifile, MPI_IO_DATA_lag_bubbles, 21*tot_data, & - mpi_p, status, ierr) - do i = 1, tot_data - id = int(MPI_IO_DATA_lag_bubbles(i, 1)) - inputvals(1:20) = MPI_IO_DATA_lag_bubbles(i, 2:21) - indomain = particle_in_domain(inputvals(1:3)) - if (indomain .and. (id > 0)) then - bub_id = bub_id + 1 - nBubs = bub_id ! local number of bubbles - lag_id(bub_id, 1) = id ! global ID - lag_id(bub_id, 2) = bub_id ! local ID - mtn_pos(bub_id, 1:3, 1) = inputvals(1:3) - mtn_posPrev(bub_id, 1:3, 1) = inputvals(4:6) - mtn_vel(bub_id, 1:3, 1) = inputvals(7:9) - intfc_rad(bub_id, 1) = inputvals(10) - intfc_vel(bub_id, 1) = inputvals(11) - bub_R0(bub_id) = inputvals(12) - Rmax_stats(bub_id) = inputvals(13) - Rmin_stats(bub_id) = inputvals(14) - bub_dphidt(bub_id) = inputvals(15) - gas_p(bub_id, 1) = inputvals(16) - gas_mv(bub_id, 1) = inputvals(17) - gas_mg(bub_id) = inputvals(18) - gas_betaT(bub_id) = inputvals(19) - gas_betaC(bub_id) = inputvals(20) - cell = -buff_size - call s_locate_cell(mtn_pos(bub_id, 1:3, 1), cell, mtn_s(bub_id, 1:3, 1)) - end if + + ! Skip extended header + disp = int(sizeof(file_tot_part) + 2*sizeof(file_time) + sizeof(file_num_procs) + & + file_num_procs*sizeof(proc_bubble_counts(1)), MPI_OFFSET_KIND) + call MPI_FILE_SET_VIEW(ifile, disp, mpi_p, view, 'native', mpi_info_int, ierr) + + call MPI_FILE_READ_ALL(ifile, MPI_IO_DATA_lag_bubbles, & + lag_io_vars*bub_id, mpi_p, status, ierr) + + call MPI_FILE_CLOSE(ifile, ierr) + call MPI_TYPE_FREE(view, ierr) + + nBubs = bub_id + + do i = 1, bub_id + lag_id(i, 1) = int(MPI_IO_DATA_lag_bubbles(i, 1)) + mtn_pos(i, 1:3, 1) = MPI_IO_DATA_lag_bubbles(i, 2:4) + mtn_posPrev(i, 1:3, 1) = MPI_IO_DATA_lag_bubbles(i, 5:7) + mtn_vel(i, 1:3, 1) = MPI_IO_DATA_lag_bubbles(i, 8:10) + intfc_rad(i, 1) = MPI_IO_DATA_lag_bubbles(i, 11) + intfc_vel(i, 1) = MPI_IO_DATA_lag_bubbles(i, 12) + bub_R0(i) = MPI_IO_DATA_lag_bubbles(i, 13) + Rmax_stats(i) = MPI_IO_DATA_lag_bubbles(i, 14) + Rmin_stats(i) = MPI_IO_DATA_lag_bubbles(i, 15) + bub_dphidt(i) = MPI_IO_DATA_lag_bubbles(i, 16) + gas_p(i, 1) = MPI_IO_DATA_lag_bubbles(i, 17) + gas_mv(i, 1) = MPI_IO_DATA_lag_bubbles(i, 18) + gas_mg(i) = MPI_IO_DATA_lag_bubbles(i, 19) + gas_betaT(i) = MPI_IO_DATA_lag_bubbles(i, 20) + gas_betaC(i) = MPI_IO_DATA_lag_bubbles(i, 21) + cell = -buff_size + call s_locate_cell(mtn_pos(i, 1:3, 1), cell, mtn_s(i, 1:3, 1)) end do + deallocate (MPI_IO_DATA_lag_bubbles) + + else + nBubs = 0 + + call MPI_TYPE_CONTIGUOUS(0, mpi_p, view, ierr) + call MPI_TYPE_COMMIT(view, ierr) + + call MPI_FILE_OPEN(MPI_COMM_WORLD, file_loc, MPI_MODE_RDONLY, & + mpi_info_int, ifile, ierr) + + ! Skip extended header + disp = int(sizeof(file_tot_part) + 2*sizeof(file_time) + sizeof(file_num_procs) + & + file_num_procs*sizeof(proc_bubble_counts(1)), MPI_OFFSET_KIND) + call MPI_FILE_SET_VIEW(ifile, disp, mpi_p, view, 'native', mpi_info_int, ierr) + + call MPI_FILE_READ_ALL(ifile, dummy, 0, mpi_p, status, ierr) + + call MPI_FILE_CLOSE(ifile, ierr) + call MPI_TYPE_FREE(view, ierr) end if - call MPI_FILE_CLOSE(ifile, ierr) + + if (proc_rank == 0) then + write (*, '(A,I0,A,I0)') 'Read ', file_tot_part, ' particles from restart file at t_step = ', save_count + write (*, '(A,E15.7,A,E15.7)') 'Restart time = ', mytime, ', dt = ', dt + end if + + deallocate (proc_bubble_counts) #endif end subroutine s_restart_bubbles @@ -514,15 +622,15 @@ contains real(wp) :: myR_m, mygamma_m, myPb, myMass_n, myMass_v real(wp) :: myR, myV, myBeta_c, myBeta_t, myR0, myPbdot, myMvdot real(wp) :: myPinf, aux1, aux2, myCson, myRho - real(wp) :: gamma, pi_inf, qv + real(wp), dimension(3) :: myPos, myVel + real(wp) :: gamma, pi_inf, qv, f_b real(wp), dimension(contxe) :: myalpha_rho, myalpha real(wp), dimension(2) :: Re integer, dimension(3) :: cell - integer :: adap_dt_stop_max, adap_dt_stop !< Fail-safe exit if max iteration count reached real(wp) :: dmalf, dmntait, dmBtait, dm_bub_adv_src, dm_divu !< Dummy variables for unified subgrid bubble subroutines - integer :: i, k, l + integer :: i, j, k, l call nvtxStartRange("LAGRANGE-BUBBLE-DYNAMICS") @@ -550,11 +658,10 @@ contains ! Radial motion model adap_dt_stop_max = 0 - $:GPU_PARALLEL_LOOP(private='[k,myalpha_rho,myalpha,Re,cell]', & + $:GPU_PARALLEL_LOOP(private='[k,myalpha_rho,myalpha,Re,cell,myPos,myVel]', & & reduction='[[adap_dt_stop_max]]',reductionOp='[MAX]', & & copy='[adap_dt_stop_max]',copyin='[stage]') do k = 1, nBubs - ! Keller-Miksis model ! Current bubble state myPb = gas_p(k, 2) @@ -565,6 +672,8 @@ contains myBeta_c = gas_betaC(k) myBeta_t = gas_betaT(k) myR0 = bub_R0(k) + myPos = mtn_pos(k, :, 2) + myVel = mtn_vel(k, :, 2) ! Vapor and heat fluxes call s_vflux(myR, myV, myPb, myMass_v, k, myVapFlux, myMass_n, myBeta_c, myR_m, mygamma_m) @@ -588,11 +697,11 @@ contains adap_dt_stop = 0 if (adap_dt) then - call s_advance_step(myRho, myPinf, myR, myV, myR0, myPb, myPbdot, dmalf, & dmntait, dmBtait, dm_bub_adv_src, dm_divu, & k, myMass_v, myMass_n, myBeta_c, & - myBeta_t, myCson, adap_dt_stop) + myBeta_t, myCson, adap_dt_stop, Re(1), & + myPos, myVel, cell, q_prim_vf) ! Update bubble state intfc_rad(k, 1) = myR @@ -600,8 +709,17 @@ contains gas_p(k, 1) = myPb gas_mv(k, 1) = myMass_v - else + if (lag_params%vel_model == 1) then + mtn_posPrev(k, :, 1) = mtn_pos(k, :, 2) + mtn_pos(k, :, 1) = myPos + mtn_vel(k, :, 1) = myVel + elseif (lag_params%vel_model == 2) then + mtn_posPrev(k, :, 1) = mtn_pos(k, :, 2) + mtn_pos(k, :, 1) = myPos + mtn_vel(k, :, 1) = myVel + end if + else ! Radial acceleration from bubble models intfc_dveldt(k, stage) = f_rddot(myRho, myPinf, myR, myV, myR0, & myPb, myPbdot, dmalf, dmntait, dmBtait, & @@ -611,6 +729,30 @@ contains gas_dmvdt(k, stage) = myMvdot gas_dpdt(k, stage) = myPbdot + do l = 1, num_dims + if (lag_params%vel_model == 1) then + mtn_dposdt(k, l, stage) = f_interpolate_velocity(myPos(l), & + cell, l, q_prim_vf) + mtn_dveldt(k, l, stage) = 0._wp + elseif (lag_params%vel_model == 2) then + mtn_dposdt(k, l, stage) = myVel(l) + f_b = f_get_bubble_force(myPos(l), & + myR, myV, myVel, & + myMass_n, myMass_v, & + Re(1), myRho, cell, l, q_prim_vf) + mtn_dveldt(k, l, stage) = f_b / (myMass_n + myMass_v) + elseif (lag_params%vel_model == 3) then + mtn_dposdt(k, l, stage) = myVel(l) + f_b = f_get_bubble_force(myPos(l), & + myR, myV, myVel, & + myMass_n, myMass_v, & + Re(1), myRho, cell, l, q_prim_vf) + mtn_dveldt(k, l, stage) = 2._wp * f_b / (myMass_n + myMass_v) - 3._wp * myV * myVel(l) / myR + else + mtn_dposdt(k, l, stage) = 0._wp + mtn_dveldt(k, l, stage) = 0._wp + end if + end do end if adap_dt_stop_max = max(adap_dt_stop_max, adap_dt_stop) @@ -619,14 +761,9 @@ contains if (adap_dt .and. adap_dt_stop_max > 0) call s_mpi_abort("Adaptive time stepping failed to converge.") - ! Bubbles remain in a fixed position - $:GPU_PARALLEL_LOOP(collapse=2, private='[k]', copyin='[stage]') - do k = 1, nBubs - do l = 1, 3 - mtn_dposdt(k, l, stage) = 0._wp - mtn_dveldt(k, l, stage) = 0._wp - end do - end do + if (adap_dt .and. lag_params%vel_model > 0) then + call s_enforce_EL_bubbles_boundary_conditions(q_prim_vf, dest=1) + end if call nvtxEndRange @@ -649,6 +786,7 @@ contains if (lag_params%solver_approach == 2) then + ! (q / (1 - beta)) * d(beta)/dt source if (p == 0) then $:GPU_PARALLEL_LOOP(collapse=4) do k = 0, p @@ -686,6 +824,7 @@ contains call s_gradient_dir(q_prim_vf(E_idx), q_beta%vf(3), l) + ! (beta / (1 - beta)) * dP/dl source $:GPU_PARALLEL_LOOP(collapse=3) do k = 0, p do j = 0, n @@ -712,6 +851,7 @@ contains call s_gradient_dir(q_beta%vf(3), q_beta%vf(4), l) + ! (beta / (1 - beta)) * d(Pu)/dl source $:GPU_PARALLEL_LOOP(collapse=3) do k = 0, p do j = 0, n @@ -725,7 +865,6 @@ contains end do end do end do - end if end subroutine s_compute_bubbles_EL_source @@ -807,7 +946,7 @@ contains !! @param f_pinfl Driving pressure !! @param cell Bubble cell !! @param Romega Control volume radius - pure subroutine s_get_pinf(bub_id, q_prim_vf, ptype, f_pinfl, cell, preterm1, term2, Romega) + impure subroutine s_get_pinf(bub_id, q_prim_vf, ptype, f_pinfl, cell, preterm1, term2, Romega) $:GPU_ROUTINE(function_name='s_get_pinf',parallelism='[seq]', & & cray_inline=True) @@ -817,7 +956,8 @@ contains integer, dimension(3), intent(out) :: cell real(wp), intent(out), optional :: preterm1, term2, Romega - real(wp), dimension(3) :: scoord, psi + real(wp), dimension(3) :: scoord, psi_pos, psi_x, psi_y, psi_z + real(wp) :: xi, eta, zeta real(wp) :: dc, vol, aux real(wp) :: volgas, term1, Rbeq, denom real(wp) :: charvol, charpres, charvol2, charpres2 @@ -826,83 +966,167 @@ contains integer :: smearGrid, smearGridz logical :: celloutside - scoord = mtn_s(bub_id, 1:3, 2) f_pinfl = 0._wp - !< Find current bubble cell - cell(:) = int(scoord(:)) - $:GPU_LOOP(parallelism='[seq]') - do i = 1, num_dims - if (scoord(i) < 0._wp) cell(i) = cell(i) - 1 - end do + if (lag_params%vel_model > 0) then + cell = fd_number - buff_size + call s_locate_cell(mtn_pos(bub_id, 1:3, 2), cell, mtn_s(bub_id, 1:3, 2)) + scoord = mtn_s(bub_id, 1:3, 2) + else + scoord = mtn_s(bub_id, 1:3, 2) + cell(:) = int(scoord(:)) + $:GPU_LOOP(parallelism='[seq]') + do i = 1, num_dims + if (scoord(i) < 0._wp) cell(i) = cell(i) - 1 + end do + end if if ((lag_params%cluster_type == 1)) then !< Getting p_cell in terms of only the current cell by interpolation - !< Getting the cell volulme as Omega - if (p > 0) then - vol = dx(cell(1))*dy(cell(2))*dz(cell(3)) - else - if (cyl_coord) then - vol = dx(cell(1))*dy(cell(2))*y_cc(cell(2))*2._wp*pi + if (fd_order == 2) then ! Bilinear interpolation + + if (p > 0) then + vol = dx(cell(1))*dy(cell(2))*dz(cell(3)) else - vol = dx(cell(1))*dy(cell(2))*lag_params%charwidth + if (cyl_coord) then + vol = dx(cell(1))*dy(cell(2))*y_cc(cell(2))*2._wp*pi + else + vol = dx(cell(1))*dy(cell(2))*lag_params%charwidth + end if end if - end if - !< Obtain bilinear interpolation coefficients, based on the current location of the bubble. - psi(1) = (scoord(1) - real(cell(1)))*dx(cell(1)) + x_cb(cell(1) - 1) - if (cell(1) == (m + buff_size)) then - cell(1) = cell(1) - 1 - psi(1) = 1._wp - else if (cell(1) == (-buff_size)) then - psi(1) = 0._wp - else - if (psi(1) < x_cc(cell(1))) cell(1) = cell(1) - 1 - psi(1) = abs((psi(1) - x_cc(cell(1)))/(x_cc(cell(1) + 1) - x_cc(cell(1)))) - end if + !< Obtain bilinear interpolation coefficients, based on the current location of the bubble. + psi_pos(1) = (scoord(1) - real(cell(1)))*dx(cell(1)) + x_cb(cell(1) - 1) + psi_pos(1) = abs((psi_pos(1) - x_cc(cell(1)))/(x_cc(cell(1) + 1) - x_cc(cell(1)))) - psi(2) = (scoord(2) - real(cell(2)))*dy(cell(2)) + y_cb(cell(2) - 1) - if (cell(2) == (n + buff_size)) then - cell(2) = cell(2) - 1 - psi(2) = 1._wp - else if (cell(2) == (-buff_size)) then - psi(2) = 0._wp - else - if (psi(2) < y_cc(cell(2))) cell(2) = cell(2) - 1 - psi(2) = abs((psi(2) - y_cc(cell(2)))/(y_cc(cell(2) + 1) - y_cc(cell(2)))) - end if + psi_pos(2) = (scoord(2) - real(cell(2)))*dy(cell(2)) + y_cb(cell(2) - 1) + psi_pos(2) = abs((psi_pos(2) - y_cc(cell(2)))/(y_cc(cell(2) + 1) - y_cc(cell(2)))) - if (p > 0) then - psi(3) = (scoord(3) - real(cell(3)))*dz(cell(3)) + z_cb(cell(3) - 1) - if (cell(3) == (p + buff_size)) then - cell(3) = cell(3) - 1 - psi(3) = 1._wp - else if (cell(3) == (-buff_size)) then - psi(3) = 0._wp + if (p > 0) then + psi_pos(3) = (scoord(3) - real(cell(3)))*dz(cell(3)) + z_cb(cell(3) - 1) + psi_pos(3) = abs((psi_pos(3) - z_cc(cell(3)))/(z_cc(cell(3) + 1) - z_cc(cell(3)))) else - if (psi(3) < z_cc(cell(3))) cell(3) = cell(3) - 1 - psi(3) = abs((psi(3) - z_cc(cell(3)))/(z_cc(cell(3) + 1) - z_cc(cell(3)))) + psi_pos(3) = 0._wp end if - else - psi(3) = 0._wp - end if - !< Perform bilinear interpolation - if (p == 0) then !2D - f_pinfl = q_prim_vf(E_idx)%sf(cell(1), cell(2), cell(3))*(1._wp - psi(1))*(1._wp - psi(2)) - f_pinfl = f_pinfl + q_prim_vf(E_idx)%sf(cell(1) + 1, cell(2), cell(3))*psi(1)*(1._wp - psi(2)) - f_pinfl = f_pinfl + q_prim_vf(E_idx)%sf(cell(1) + 1, cell(2) + 1, cell(3))*psi(1)*psi(2) - f_pinfl = f_pinfl + q_prim_vf(E_idx)%sf(cell(1), cell(2) + 1, cell(3))*(1._wp - psi(1))*psi(2) - else !3D - f_pinfl = q_prim_vf(E_idx)%sf(cell(1), cell(2), cell(3))*(1._wp - psi(1))*(1._wp - psi(2))*(1._wp - psi(3)) - f_pinfl = f_pinfl + q_prim_vf(E_idx)%sf(cell(1) + 1, cell(2), cell(3))*psi(1)*(1._wp - psi(2))*(1._wp - psi(3)) - f_pinfl = f_pinfl + q_prim_vf(E_idx)%sf(cell(1) + 1, cell(2) + 1, cell(3))*psi(1)*psi(2)*(1._wp - psi(3)) - f_pinfl = f_pinfl + q_prim_vf(E_idx)%sf(cell(1), cell(2) + 1, cell(3))*(1._wp - psi(1))*psi(2)*(1._wp - psi(3)) - f_pinfl = f_pinfl + q_prim_vf(E_idx)%sf(cell(1), cell(2), cell(3) + 1)*(1._wp - psi(1))*(1._wp - psi(2))*psi(3) - f_pinfl = f_pinfl + q_prim_vf(E_idx)%sf(cell(1) + 1, cell(2), cell(3) + 1)*psi(1)*(1._wp - psi(2))*psi(3) - f_pinfl = f_pinfl + q_prim_vf(E_idx)%sf(cell(1) + 1, cell(2) + 1, cell(3) + 1)*psi(1)*psi(2)*psi(3) - f_pinfl = f_pinfl + q_prim_vf(E_idx)%sf(cell(1), cell(2) + 1, cell(3) + 1)*(1._wp - psi(1))*psi(2)*psi(3) + ! Calculate bilinear basis functions for each direction + ! For normalized coordinate xi in [0, 1], the two basis functions are: + ! phi_0(xi) = 1 - xi, phi_1(xi) = xi + + ! X-direction basis functions + psi_x(1) = 1._wp - psi_pos(1) ! Left basis function + psi_x(2) = psi_pos(1) ! Right basis function + + ! Y-direction basis functions + psi_y(1) = 1._wp - psi_pos(2) ! Left basis function + psi_y(2) = psi_pos(2) ! Right basis function + + if (p > 0) then + ! Z-direction basis functions + psi_z(1) = 1._wp - psi_pos(3) ! Left basis function + psi_z(2) = psi_pos(3) ! Right basis function + else + psi_z(1) = 1._wp + psi_z(2) = 0._wp + end if + + !< Perform bilinear interpolation + f_pinfl = 0._wp + + if (p == 0) then !2D - 4 point interpolation (2x2) + do j = 1, 2 + do i = 1, 2 + f_pinfl = f_pinfl + q_prim_vf(E_idx)%sf(cell(1) + i - 1, cell(2) + j - 1, cell(3))* & + psi_x(i)*psi_y(j) + end do + end do + else !3D - 8 point interpolation (2x2x2) + do k = 1, 2 + do j = 1, 2 + do i = 1, 2 + f_pinfl = f_pinfl + q_prim_vf(E_idx)%sf(cell(1) + i - 1, cell(2) + j - 1, cell(3) + k - 1)* & + psi_x(i)*psi_y(j)*psi_z(k) + end do + end do + end do + end if + + elseif (fd_order == 4) then ! Biquadratic interpolation + + if (p > 0) then + vol = dx(cell(1))*dy(cell(2))*dz(cell(3)) + else + if (cyl_coord) then + vol = dx(cell(1))*dy(cell(2))*y_cc(cell(2))*2._wp*pi + else + vol = dx(cell(1))*dy(cell(2))*lag_params%charwidth + end if + end if + + !< Obtain biquadratic interpolation coefficients, based on the current location of the bubble. + ! For biquadratic interpolation, we need coefficients for 3 points in each direction + psi_pos(1) = (scoord(1) - real(cell(1)))*dx(cell(1)) + x_cb(cell(1) - 1) + psi_pos(1) = (psi_pos(1) - x_cc(cell(1)))/(x_cc(cell(1) + 1) - x_cc(cell(1))) + + psi_pos(2) = (scoord(2) - real(cell(2)))*dy(cell(2)) + y_cb(cell(2) - 1) + psi_pos(2) = (psi_pos(2) - y_cc(cell(2)))/(y_cc(cell(2) + 1) - y_cc(cell(2))) + + if (p > 0) then + psi_pos(3) = (scoord(3) - real(cell(3)))*dz(cell(3)) + z_cb(cell(3) - 1) + psi_pos(3) = (psi_pos(3) - z_cc(cell(3)))/(z_cc(cell(3) + 1) - z_cc(cell(3))) + else + psi_pos(3) = 0._wp + end if + + ! Calculate biquadratic basis functions for each direction + ! For normalized coordinate xi in [-1, 1], the three basis functions are: + ! phi_0(xi) = xi*(xi-1)/2, phi_1(xi) = (1-xi)*(1+xi), phi_2(xi) = xi*(xi+1)/2 + + ! X-direction basis functions + xi = 2._wp*psi_pos(1) - 1._wp ! Convert to [-1, 1] range + psi_x(1) = xi*(xi - 1._wp)/2._wp ! Left basis function + psi_x(2) = (1._wp - xi)*(1._wp + xi) ! Center basis function + psi_x(3) = xi*(xi + 1._wp)/2._wp ! Right basis function + + ! Y-direction basis functions + eta = 2._wp*psi_pos(2) - 1._wp ! Convert to [-1, 1] range + psi_y(1) = eta*(eta - 1._wp)/2._wp ! Left basis function + psi_y(2) = (1._wp - eta)*(1._wp + eta) ! Center basis function + psi_y(3) = eta*(eta + 1._wp)/2._wp ! Right basis function + + if (p > 0) then + ! Z-direction basis functions + zeta = 2._wp*psi_pos(3) - 1._wp ! Convert to [-1, 1] range + psi_z(1) = zeta*(zeta - 1._wp)/2._wp ! Left basis function + psi_z(2) = (1._wp - zeta)*(1._wp + zeta) ! Center basis function + psi_z(3) = zeta*(zeta + 1._wp)/2._wp ! Right basis function + else + psi_z(1) = 0._wp + psi_z(2) = 1._wp + psi_z(3) = 0._wp + end if + + !< Perform biquadratic interpolation + f_pinfl = 0._wp + + if (p == 0) then !2D - 9 point interpolation (3x3) + do j = 1, 3 + do i = 1, 3 + f_pinfl = f_pinfl + q_prim_vf(E_idx)%sf(cell(1) + i - 2, cell(2) + j - 2, cell(3))* & + psi_x(i)*psi_y(j) + end do + end do + else !3D - 27 point interpolation (3x3x3) + do k = 1, 3 + do j = 1, 3 + do i = 1, 3 + f_pinfl = f_pinfl + q_prim_vf(E_idx)%sf(cell(1) + i - 2, cell(2) + j - 2, cell(3) + k - 2)* & + psi_x(i)*psi_y(j)*psi_z(k) + end do + end do + end do + end if end if !R_Omega @@ -936,26 +1160,32 @@ contains !< check if the current cell is outside the computational domain or not (including ghost cells) celloutside = .false. if (num_dims == 2) then - if ((cellaux(1) < -buff_size) .or. (cellaux(2) < -buff_size)) then + if ((cellaux(1) < fd_number - buff_size) .or. & + (cellaux(2) < fd_number - buff_size)) then celloutside = .true. end if - if (cyl_coord .and. y_cc(cellaux(2)) < 0._wp) then + if (cyl_coord .and. cellaux(2) < 0) then celloutside = .true. end if - if ((cellaux(2) > n + buff_size) .or. (cellaux(1) > m + buff_size)) then + if ((cellaux(2) > n + buff_size - fd_number) .or. & + (cellaux(1) > m + buff_size - fd_number)) then celloutside = .true. end if else - if ((cellaux(3) < -buff_size) .or. (cellaux(1) < -buff_size) .or. (cellaux(2) < -buff_size)) then + if ((cellaux(3) < fd_number - buff_size) .or. & + (cellaux(1) < fd_number - buff_size) .or. & + (cellaux(2) < fd_number - buff_size)) then celloutside = .true. end if - if ((cellaux(3) > p + buff_size) .or. (cellaux(2) > n + buff_size) .or. (cellaux(1) > m + buff_size)) then + if ((cellaux(3) > p + buff_size - fd_number) .or. & + (cellaux(2) > n + buff_size - fd_number) .or. & + (cellaux(1) > m + buff_size - fd_number)) then celloutside = .true. end if end if if (.not. celloutside) then - if (cyl_coord .and. (p == 0) .and. (y_cc(cellaux(2)) < 0._wp)) then + if (cyl_coord .and. (p == 0) .and. (cellaux(2) < 0)) then celloutside = .true. end if end if @@ -1018,8 +1248,9 @@ contains !> This subroutine updates the Lagrange variables using the tvd RK time steppers. !! The time derivative of the bubble variables must be stored at every stage to avoid precision errors. !! @param stage Current tvd RK stage - impure subroutine s_update_lagrange_tdv_rk(stage) + impure subroutine s_update_lagrange_tdv_rk(q_prim_vf, stage) + type(scalar_field), dimension(sys_size), intent(in) :: q_prim_vf integer, intent(in) :: stage integer :: k @@ -1030,12 +1261,15 @@ contains !u{1} = u{n} + dt * RHS{n} intfc_rad(k, 1) = intfc_rad(k, 1) + dt*intfc_draddt(k, 1) intfc_vel(k, 1) = intfc_vel(k, 1) + dt*intfc_dveldt(k, 1) + mtn_posPrev(k, 1:3, 1) = mtn_pos(k, 1:3, 1) mtn_pos(k, 1:3, 1) = mtn_pos(k, 1:3, 1) + dt*mtn_dposdt(k, 1:3, 1) mtn_vel(k, 1:3, 1) = mtn_vel(k, 1:3, 1) + dt*mtn_dveldt(k, 1:3, 1) gas_p(k, 1) = gas_p(k, 1) + dt*gas_dpdt(k, 1) gas_mv(k, 1) = gas_mv(k, 1) + dt*gas_dmvdt(k, 1) end do + if (lag_params%vel_model > 0) call s_enforce_EL_bubbles_boundary_conditions(q_prim_vf, dest=1) + call s_transfer_data_to_tmp() call s_write_void_evol(mytime) if (lag_params%write_bubbles_stats) call s_calculate_lag_bubble_stats() @@ -1052,28 +1286,33 @@ contains !u{1} = u{n} + dt * RHS{n} intfc_rad(k, 2) = intfc_rad(k, 1) + dt*intfc_draddt(k, 1) intfc_vel(k, 2) = intfc_vel(k, 1) + dt*intfc_dveldt(k, 1) + mtn_posPrev(k, 1:3, 2) = mtn_pos(k, 1:3, 1) mtn_pos(k, 1:3, 2) = mtn_pos(k, 1:3, 1) + dt*mtn_dposdt(k, 1:3, 1) mtn_vel(k, 1:3, 2) = mtn_vel(k, 1:3, 1) + dt*mtn_dveldt(k, 1:3, 1) gas_p(k, 2) = gas_p(k, 1) + dt*gas_dpdt(k, 1) gas_mv(k, 2) = gas_mv(k, 1) + dt*gas_dmvdt(k, 1) end do + if (lag_params%vel_model > 0) call s_enforce_EL_bubbles_boundary_conditions(q_prim_vf, dest=2) + elseif (stage == 2) then $:GPU_PARALLEL_LOOP(private='[k]') do k = 1, nBubs !u{1} = u{n} + (1/2) * dt * (RHS{n} + RHS{1}) intfc_rad(k, 1) = intfc_rad(k, 1) + dt*(intfc_draddt(k, 1) + intfc_draddt(k, 2))/2._wp intfc_vel(k, 1) = intfc_vel(k, 1) + dt*(intfc_dveldt(k, 1) + intfc_dveldt(k, 2))/2._wp + mtn_posPrev(k, 1:3, 1) = mtn_pos(k, 1:3, 2) mtn_pos(k, 1:3, 1) = mtn_pos(k, 1:3, 1) + dt*(mtn_dposdt(k, 1:3, 1) + mtn_dposdt(k, 1:3, 2))/2._wp mtn_vel(k, 1:3, 1) = mtn_vel(k, 1:3, 1) + dt*(mtn_dveldt(k, 1:3, 1) + mtn_dveldt(k, 1:3, 2))/2._wp gas_p(k, 1) = gas_p(k, 1) + dt*(gas_dpdt(k, 1) + gas_dpdt(k, 2))/2._wp gas_mv(k, 1) = gas_mv(k, 1) + dt*(gas_dmvdt(k, 1) + gas_dmvdt(k, 2))/2._wp end do + if (lag_params%vel_model > 0) call s_enforce_EL_bubbles_boundary_conditions(q_prim_vf, dest=1) + call s_transfer_data_to_tmp() call s_write_void_evol(mytime) if (lag_params%write_bubbles_stats) call s_calculate_lag_bubble_stats() - if (lag_params%write_bubbles) then $:GPU_UPDATE(host='[gas_p,gas_mv,intfc_rad,intfc_vel]') call s_write_lag_particles(mytime) @@ -1088,35 +1327,45 @@ contains !u{1} = u{n} + dt * RHS{n} intfc_rad(k, 2) = intfc_rad(k, 1) + dt*intfc_draddt(k, 1) intfc_vel(k, 2) = intfc_vel(k, 1) + dt*intfc_dveldt(k, 1) + mtn_posPrev(k, 1:3, 2) = mtn_pos(k, 1:3, 1) mtn_pos(k, 1:3, 2) = mtn_pos(k, 1:3, 1) + dt*mtn_dposdt(k, 1:3, 1) mtn_vel(k, 1:3, 2) = mtn_vel(k, 1:3, 1) + dt*mtn_dveldt(k, 1:3, 1) gas_p(k, 2) = gas_p(k, 1) + dt*gas_dpdt(k, 1) gas_mv(k, 2) = gas_mv(k, 1) + dt*gas_dmvdt(k, 1) end do + if (lag_params%vel_model > 0) call s_enforce_EL_bubbles_boundary_conditions(q_prim_vf, dest=2) + elseif (stage == 2) then $:GPU_PARALLEL_LOOP(private='[k]') do k = 1, nBubs !u{2} = u{n} + (1/4) * dt * [RHS{n} + RHS{1}] intfc_rad(k, 2) = intfc_rad(k, 1) + dt*(intfc_draddt(k, 1) + intfc_draddt(k, 2))/4._wp intfc_vel(k, 2) = intfc_vel(k, 1) + dt*(intfc_dveldt(k, 1) + intfc_dveldt(k, 2))/4._wp + mtn_posPrev(k, 1:3, 2) = mtn_pos(k, 1:3, 2) mtn_pos(k, 1:3, 2) = mtn_pos(k, 1:3, 1) + dt*(mtn_dposdt(k, 1:3, 1) + mtn_dposdt(k, 1:3, 2))/4._wp mtn_vel(k, 1:3, 2) = mtn_vel(k, 1:3, 1) + dt*(mtn_dveldt(k, 1:3, 1) + mtn_dveldt(k, 1:3, 2))/4._wp gas_p(k, 2) = gas_p(k, 1) + dt*(gas_dpdt(k, 1) + gas_dpdt(k, 2))/4._wp gas_mv(k, 2) = gas_mv(k, 1) + dt*(gas_dmvdt(k, 1) + gas_dmvdt(k, 2))/4._wp end do + + if (lag_params%vel_model > 0) call s_enforce_EL_bubbles_boundary_conditions(q_prim_vf, dest=2) + elseif (stage == 3) then $:GPU_PARALLEL_LOOP(private='[k]') do k = 1, nBubs !u{n+1} = u{n} + (2/3) * dt * [(1/4)* RHS{n} + (1/4)* RHS{1} + RHS{2}] intfc_rad(k, 1) = intfc_rad(k, 1) + (2._wp/3._wp)*dt*(intfc_draddt(k, 1)/4._wp + intfc_draddt(k, 2)/4._wp + intfc_draddt(k, 3)) intfc_vel(k, 1) = intfc_vel(k, 1) + (2._wp/3._wp)*dt*(intfc_dveldt(k, 1)/4._wp + intfc_dveldt(k, 2)/4._wp + intfc_dveldt(k, 3)) + mtn_posPrev(k, 1:3, 1) = mtn_pos(k, 1:3, 2) mtn_pos(k, 1:3, 1) = mtn_pos(k, 1:3, 1) + (2._wp/3._wp)*dt*(mtn_dposdt(k, 1:3, 1)/4._wp + mtn_dposdt(k, 1:3, 2)/4._wp + mtn_dposdt(k, 1:3, 3)) mtn_vel(k, 1:3, 1) = mtn_vel(k, 1:3, 1) + (2._wp/3._wp)*dt*(mtn_dveldt(k, 1:3, 1)/4._wp + mtn_dveldt(k, 1:3, 2)/4._wp + mtn_dveldt(k, 1:3, 3)) gas_p(k, 1) = gas_p(k, 1) + (2._wp/3._wp)*dt*(gas_dpdt(k, 1)/4._wp + gas_dpdt(k, 2)/4._wp + gas_dpdt(k, 3)) gas_mv(k, 1) = gas_mv(k, 1) + (2._wp/3._wp)*dt*(gas_dmvdt(k, 1)/4._wp + gas_dmvdt(k, 2)/4._wp + gas_dmvdt(k, 3)) end do + if (lag_params%vel_model > 0) call s_enforce_EL_bubbles_boundary_conditions(q_prim_vf, dest=1) + call s_transfer_data_to_tmp() call s_write_void_evol(mytime) if (lag_params%write_bubbles_stats) call s_calculate_lag_bubble_stats() @@ -1127,16 +1376,221 @@ contains end if end if - end if end subroutine s_update_lagrange_tdv_rk + !> This subroutine enforces reflective and wall boundary conditions for EL bubbles + !! @param dest Destination for the bubble position update + impure subroutine s_enforce_EL_bubbles_boundary_conditions(q_prim_vf, dest) + + type(scalar_field), dimension(sys_size), intent(in) :: q_prim_vf + integer, intent(in) :: dest + integer :: k, i, q + integer :: patch_id, newBubs, new_idx + real(wp) :: offset + integer, dimension(3) :: cell + + $:GPU_PARALLEL_LOOP(private='[cell]') + do k = 1, nBubs + keep_bubble(k) = 1 + wrap_bubble_loc(k,:) = 0 + wrap_bubble_dir(k,:) = 0 + + ! Relocate bubbles at solid boundaries and delete bubbles that leave + ! buffer regions + if (any(bc_x%beg == (/BC_REFLECTIVE, BC_CHAR_SLIP_WALL, BC_SLIP_WALL, BC_NO_SLIP_WALL/)) & + .and. mtn_pos(k, 1, dest) < x_cb(-1) + intfc_rad(k, dest)) then + mtn_pos(k, 1, dest) = x_cb(-1) + intfc_rad(k, dest) + elseif (any(bc_x%end == (/BC_REFLECTIVE, BC_CHAR_SLIP_WALL, BC_SLIP_WALL, BC_NO_SLIP_WALL/)) & + .and. mtn_pos(k, 1, dest) > x_cb(m) - intfc_rad(k, dest)) then + mtn_pos(k, 1, dest) = x_cb(m) - intfc_rad(k, dest) + elseif (bc_x%beg == BC_PERIODIC .and. mtn_pos(k, 1, dest) < pcomm_coords(1)%beg .and. & + mtn_posPrev(k, 1, dest) > pcomm_coords(1)%beg) then + wrap_bubble_dir(k,1) = 1 + wrap_bubble_loc(k,1) = -1 + elseif (bc_x%end == BC_PERIODIC .and. mtn_pos(k, 1, dest) > pcomm_coords(1)%end .and. & + mtn_posPrev(k, 1, dest) < pcomm_coords(1)%end) then + wrap_bubble_dir(k,1) = 1 + wrap_bubble_loc(k,1) = 1 + elseif (mtn_pos(k, 1, dest) >= x_cb(m + buff_size - fd_number)) then + keep_bubble(k) = 0 + elseif (mtn_pos(k, 1, dest) < x_cb(fd_number - buff_size - 1)) then + keep_bubble(k) = 0 + end if + + if (any(bc_y%beg == (/BC_REFLECTIVE, BC_CHAR_SLIP_WALL, BC_SLIP_WALL, BC_NO_SLIP_WALL/)) & + .and. mtn_pos(k, 2, dest) < y_cb(-1) + intfc_rad(k, dest)) then + mtn_pos(k, 2, dest) = y_cb(-1) + intfc_rad(k, dest) + else if (any(bc_y%end == (/BC_REFLECTIVE, BC_CHAR_SLIP_WALL, BC_SLIP_WALL, BC_NO_SLIP_WALL/)) & + .and. mtn_pos(k, 2, dest) > y_cb(n) - intfc_rad(k, dest)) then + mtn_pos(k, 2, dest) = y_cb(n) - intfc_rad(k, dest) + elseif (bc_y%beg == BC_PERIODIC .and. mtn_pos(k, 2, dest) < pcomm_coords(2)%beg .and. & + mtn_posPrev(k, 2, dest) > pcomm_coords(2)%beg) then + wrap_bubble_dir(k,2) = 1 + wrap_bubble_loc(k,2) = -1 + elseif (bc_y%end == BC_PERIODIC .and. mtn_pos(k, 2, dest) > pcomm_coords(2)%end .and. & + mtn_posPrev(k, 2, dest) < pcomm_coords(2)%end) then + wrap_bubble_dir(k,2) = 1 + wrap_bubble_loc(k,2) = 1 + elseif (mtn_pos(k, 2, dest) >= y_cb(n + buff_size - fd_number)) then + keep_bubble(k) = 0 + elseif (mtn_pos(k, 2, dest) < y_cb(fd_number - buff_size - 1)) then + keep_bubble(k) = 0 + end if + + if (p > 0) then + if (any(bc_z%beg == (/BC_REFLECTIVE, BC_CHAR_SLIP_WALL, BC_SLIP_WALL, BC_NO_SLIP_WALL/)) & + .and. mtn_pos(k, 3, dest) < z_cb(-1) + intfc_rad(k, dest)) then + mtn_pos(k, 3, dest) = z_cb(-1) + intfc_rad(k, dest) + else if (any(bc_z%end == (/BC_REFLECTIVE, BC_CHAR_SLIP_WALL, BC_SLIP_WALL, BC_NO_SLIP_WALL/)) & + .and. mtn_pos(k, 3, dest) > z_cb(p) - intfc_rad(k, dest)) then + mtn_pos(k, 3, dest) = z_cb(p) - intfc_rad(k, dest) + elseif (bc_z%beg == BC_PERIODIC .and. mtn_pos(k, 3, dest) < pcomm_coords(3)%beg .and. & + mtn_posPrev(k, 3, dest) > pcomm_coords(3)%beg) then + wrap_bubble_dir(k,3) = 1 + wrap_bubble_loc(k,3) = -1 + elseif (bc_z%end == BC_PERIODIC .and. mtn_pos(k, 3, dest) > pcomm_coords(3)%end .and. & + mtn_posPrev(k, 3, dest) < pcomm_coords(3)%end) then + wrap_bubble_dir(k,3) = 1 + wrap_bubble_loc(k,3) = 1 + elseif (mtn_pos(k, 3, dest) >= z_cb(p + buff_size - fd_number)) then + keep_bubble(k) = 0 + elseif (mtn_pos(k, 3, dest) < z_cb(fd_number - buff_size - 1)) then + keep_bubble(k) = 0 + end if + end if + + if (keep_bubble(k) == 1) then + ! Remove bubbles that are no longer in a liquid + cell = fd_number - buff_size + call s_locate_cell(mtn_pos(k, 1:3, dest), cell, mtn_s(k, 1:3, dest)) + + if (q_prim_vf(advxb)%sf(cell(1), cell(2), cell(3)) < (1._wp - lag_params%valmaxvoid)) then + keep_bubble(k) = 0 + end if + + ! Move bubbles back to surface of IB + if (ib) then + cell = fd_number - buff_size + call s_locate_cell(mtn_pos(k, 1:3, dest), cell, mtn_s(k, 1:3, dest)) + + if (ib_markers%sf(cell(1), cell(2), cell(3)) /= 0) then + patch_id = ib_markers%sf(cell(1), cell(2), cell(3)) + + $:GPU_LOOP(parallelism='[seq]') + do i = 1, num_dims + mtn_pos(k, i, dest) = mtn_pos(k, i, dest) - & + levelset_norm%sf(cell(1), cell(2), cell(3), patch_id, i) & + *levelset%sf(cell(1), cell(2), cell(3), patch_id) + end do + + cell = fd_number - buff_size + call s_locate_cell(mtn_pos(k, 1:3, dest), cell, mtn_s(k, 1:3, dest)) + end if + end if + end if + end do + + call nvtxStartRange("LAG-BC") + call nvtxStartRange("LAG-BC-DEV2HOST") + $:GPU_UPDATE(host='[bub_R0, Rmax_stats, Rmin_stats, gas_mg, gas_betaT, & + & gas_betaC, bub_dphidt, lag_id, gas_p, gas_mv, intfc_rad, intfc_vel, & + & mtn_pos, mtn_posPrev, mtn_vel, mtn_s, intfc_draddt, intfc_dveldt, & + & gas_dpdt, gas_dmvdt, mtn_dposdt, mtn_dveldt, keep_bubble, nBubs, & + & wrap_bubble_dir, wrap_bubble_loc]') + call nvtxEndRange + + if (nBubs > 0) then + ! Handle deletion of bubbles leaving local domain + do k = 1, nBubs + if (k == 1) then + prefix_sum(k) = keep_bubble(k) + else + prefix_sum(k) = prefix_sum(k - 1) + keep_bubble(k) + end if + end do + + active_bubs = prefix_sum(nBubs) + + do k = 1, nBubs + if (keep_bubble(k) == 1) then + if (prefix_sum(k) /= k) then + call s_copy_lag_bubble(prefix_sum(k), k) + wrap_bubble_dir(prefix_sum(k), :) = wrap_bubble_dir(k, :) + wrap_bubble_loc(prefix_sum(k), :) = wrap_bubble_loc(k, :) + end if + end if + end do + + nBubs = active_bubs + + ! Handle periodic wrapping of bubbles on same processor + newBubs = 0 + do k = 1, nBubs + if (any(wrap_bubble_dir(k, :) == 1)) then + newBubs = newBubs + 1 + new_idx = nBubs + newBubs + call s_copy_lag_bubble(new_idx, k) + do i = 1, num_dims + if (wrap_bubble_dir(k, i) == 1) then + offset = glb_bounds(i)%end - glb_bounds(i)%beg + if (wrap_bubble_loc(k,i) == 1) then + do q = 1, 2 + mtn_pos(new_idx, i, q) = mtn_pos(new_idx, i, q) - offset + mtn_posPrev(new_idx, i, q) = mtn_posPrev(new_idx, i, q) - offset + end do + else if (wrap_bubble_loc(k,i) == -1) then + do q = 1, 2 + mtn_pos(new_idx, i, q) = mtn_pos(new_idx, i, q) + offset + mtn_posPrev(new_idx, i, q) = mtn_posPrev(new_idx, i, q) + offset + end do + end if + end if + end do + end if + end do + nBubs = nBubs + newBubs + end if + + if (run_time_info) then + call s_mpi_reduce_int_sum(nBubs, active_bubs) + if (proc_rank == 0 .and. active_bubs == 0) then + call s_mpi_abort('No bubbles remain in the domain. Simulation ending.') + end if + end if + + ! Handle MPI transfer of bubbles going to another processor's local domain + if (num_procs > 1) then + call nvtxStartRange("LAG-BC-TRANSFER-LIST") + call s_add_particles_to_transfer_list(nBubs, mtn_pos(:, :, dest), mtn_posPrev(:, :, dest)) + call nvtxEndRange + + call nvtxStartRange("LAG-BC-SENDRECV") + call s_mpi_sendrecv_particles(bub_R0, Rmax_stats, Rmin_stats, gas_mg, gas_betaT, & + gas_betaC, bub_dphidt, lag_id, gas_p, gas_mv, & + intfc_rad, intfc_vel, mtn_pos, mtn_posPrev, mtn_vel, & + mtn_s, intfc_draddt, intfc_dveldt, gas_dpdt, & + gas_dmvdt, mtn_dposdt, mtn_dveldt, lag_num_ts, nBubs, & + dest) + call nvtxEndRange + end if + + call nvtxStartRange("LAG-BC-HOST2DEV") + $:GPU_UPDATE(device='[bub_R0, Rmax_stats, Rmin_stats, gas_mg, gas_betaT, & + & gas_betaC, bub_dphidt, lag_id, gas_p, gas_mv, intfc_rad, intfc_vel, & + & mtn_pos, mtn_posPrev, mtn_vel, mtn_s, intfc_draddt, intfc_dveldt, & + & gas_dpdt, gas_dmvdt, mtn_dposdt, mtn_dveldt, nBubs]') + call nvtxEndRange + call nvtxEndRange + + end subroutine s_enforce_EL_bubbles_boundary_conditions + !> This subroutine returns the computational coordinate of the cell for the given position. !! @param pos Input coordinates !! @param cell Computational coordinate of the cell !! @param scoord Calculated particle coordinates - pure subroutine s_locate_cell(pos, cell, scoord) + impure subroutine s_locate_cell(pos, cell, scoord) real(wp), dimension(3), intent(in) :: pos real(wp), dimension(3), intent(out) :: scoord @@ -1148,7 +1602,7 @@ contains cell(1) = cell(1) - 1 end do - do while (pos(1) > x_cb(cell(1))) + do while (pos(1) >= x_cb(cell(1))) cell(1) = cell(1) + 1 end do @@ -1156,7 +1610,7 @@ contains cell(2) = cell(2) - 1 end do - do while (pos(2) > y_cb(cell(2))) + do while (pos(2) >= y_cb(cell(2))) cell(2) = cell(2) + 1 end do @@ -1164,7 +1618,7 @@ contains do while (pos(3) < z_cb(cell(3) - 1)) cell(3) = cell(3) - 1 end do - do while (pos(3) > z_cb(cell(3))) + do while (pos(3) >= z_cb(cell(3))) cell(3) = cell(3) + 1 end do end if @@ -1218,41 +1672,47 @@ contains if (p == 0 .and. cyl_coord .neqv. .true.) then ! Defining a virtual z-axis that has the same dimensions as y-axis ! defined in the input file - particle_in_domain = ((pos_part(1) < x_cb(m + buff_size)) .and. (pos_part(1) >= x_cb(-buff_size - 1)) .and. & - (pos_part(2) < y_cb(n + buff_size)) .and. (pos_part(2) >= y_cb(-buff_size - 1)) .and. & - (pos_part(3) < lag_params%charwidth/2._wp) .and. (pos_part(3) >= -lag_params%charwidth/2._wp)) + particle_in_domain = ((pos_part(1) < x_cb(m + buff_size - fd_number)) .and. & + (pos_part(1) >= x_cb(fd_number - buff_size - 1)) .and. & + (pos_part(2) < y_cb(n + buff_size - fd_number)) .and. & + (pos_part(2) >= y_cb(fd_number - buff_size - 1)) .and. & + (pos_part(3) < lag_params%charwidth/2._wp) .and. (pos_part(3) > -lag_params%charwidth/2._wp)) else ! cyl_coord - particle_in_domain = ((pos_part(1) < x_cb(m + buff_size)) .and. (pos_part(1) >= x_cb(-buff_size - 1)) .and. & - (abs(pos_part(2)) < y_cb(n + buff_size)) .and. (abs(pos_part(2)) >= max(y_cb(-buff_size - 1), 0._wp))) + particle_in_domain = ((pos_part(1) < x_cb(m + buff_size - fd_number)) .and. & + (pos_part(1) >= x_cb(fd_number - buff_size - 1)) .and. & + (abs(pos_part(2)) < y_cb(n + buff_size - fd_number)) .and. & + (abs(pos_part(2)) >= max(y_cb(fd_number - buff_size - 1), 0._wp))) end if ! 3D - if (p > 0) then - particle_in_domain = ((pos_part(1) < x_cb(m + buff_size)) .and. (pos_part(1) >= x_cb(-buff_size - 1)) .and. & - (pos_part(2) < y_cb(n + buff_size)) .and. (pos_part(2) >= y_cb(-buff_size - 1)) .and. & - (pos_part(3) < z_cb(p + buff_size)) .and. (pos_part(3) >= z_cb(-buff_size - 1))) + if (p > 1) then + particle_in_domain = ((pos_part(1) < x_cb(m + buff_size - fd_number)) .and. & + (pos_part(1) >= x_cb(fd_number - buff_size - 1)) .and. & + (pos_part(2) < y_cb(n + buff_size - fd_number)) .and. & + (pos_part(2) >= y_cb(fd_number - buff_size - 1)) .and. & + (pos_part(3) < z_cb(p + buff_size - fd_number)) .and. & + (pos_part(3) >= z_cb(fd_number - buff_size - 1))) end if - ! For symmetric boundary condition - if (bc_x%beg == BC_REFLECTIVE) then + ! For symmetric and wall boundary condition + if (any(bc_x%beg == (/BC_REFLECTIVE, BC_CHAR_SLIP_WALL, BC_SLIP_WALL, BC_NO_SLIP_WALL/))) then particle_in_domain = (particle_in_domain .and. (pos_part(1) >= x_cb(-1))) end if - if (bc_x%end == BC_REFLECTIVE) then + if (any(bc_x%end == (/BC_REFLECTIVE, BC_CHAR_SLIP_WALL, BC_SLIP_WALL, BC_NO_SLIP_WALL/))) then particle_in_domain = (particle_in_domain .and. (pos_part(1) < x_cb(m))) end if - if (bc_y%beg == BC_REFLECTIVE .and. (.not. cyl_coord)) then + if (any(bc_y%beg == (/BC_REFLECTIVE, BC_CHAR_SLIP_WALL, BC_SLIP_WALL, BC_NO_SLIP_WALL/)) .and. (.not. cyl_coord)) then particle_in_domain = (particle_in_domain .and. (pos_part(2) >= y_cb(-1))) end if - if (bc_y%end == BC_REFLECTIVE .and. (.not. cyl_coord)) then + if (any(bc_y%end == (/BC_REFLECTIVE, BC_CHAR_SLIP_WALL, BC_SLIP_WALL, BC_NO_SLIP_WALL/)) .and. (.not. cyl_coord)) then particle_in_domain = (particle_in_domain .and. (pos_part(2) < y_cb(n))) end if - if (p > 0) then - if (bc_z%beg == BC_REFLECTIVE) then + if (any(bc_z%beg == (/BC_REFLECTIVE, BC_CHAR_SLIP_WALL, BC_SLIP_WALL, BC_NO_SLIP_WALL/))) then particle_in_domain = (particle_in_domain .and. (pos_part(3) >= z_cb(-1))) end if - if (bc_z%end == BC_REFLECTIVE) then + if (any(bc_z%end == (/BC_REFLECTIVE, BC_CHAR_SLIP_WALL, BC_SLIP_WALL, BC_NO_SLIP_WALL/))) then particle_in_domain = (particle_in_domain .and. (pos_part(3) < z_cb(p))) end if end if @@ -1303,36 +1763,34 @@ contains end do end do end do - else - if (dir == 2) then - ! Gradient in y dir. - $:GPU_PARALLEL_LOOP(collapse=3) - do k = 0, p - do j = 0, n - do i = 0, m - dq%sf(i, j, k) = q%sf(i, j, k)*(dy(j + 1) - dy(j - 1)) & - + q%sf(i, j + 1, k)*(dy(j) + dy(j - 1)) & - - q%sf(i, j - 1, k)*(dy(j) + dy(j + 1)) - dq%sf(i, j, k) = dq%sf(i, j, k)/ & - ((dy(j) + dy(j - 1))*(dy(j) + dy(j + 1))) - end do + elseif (dir == 2) then + ! Gradient in y dir. + $:GPU_PARALLEL_LOOP(collapse=3) + do k = 0, p + do j = 0, n + do i = 0, m + dq%sf(i, j, k) = q%sf(i, j, k)*(dy(j + 1) - dy(j - 1)) & + + q%sf(i, j + 1, k)*(dy(j) + dy(j - 1)) & + - q%sf(i, j - 1, k)*(dy(j) + dy(j + 1)) + dq%sf(i, j, k) = dq%sf(i, j, k)/ & + ((dy(j) + dy(j - 1))*(dy(j) + dy(j + 1))) end do end do - else - ! Gradient in z dir. - $:GPU_PARALLEL_LOOP(collapse=3) - do k = 0, p - do j = 0, n - do i = 0, m - dq%sf(i, j, k) = q%sf(i, j, k)*(dz(k + 1) - dz(k - 1)) & - + q%sf(i, j, k + 1)*(dz(k) + dz(k - 1)) & - - q%sf(i, j, k - 1)*(dz(k) + dz(k + 1)) - dq%sf(i, j, k) = dq%sf(i, j, k)/ & - ((dz(k) + dz(k - 1))*(dz(k) + dz(k + 1))) - end do + end do + elseif (dir == 3) then + ! Gradient in z dir. + $:GPU_PARALLEL_LOOP(collapse=3) + do k = 0, p + do j = 0, n + do i = 0, m + dq%sf(i, j, k) = q%sf(i, j, k)*(dz(k + 1) - dz(k - 1)) & + + q%sf(i, j, k + 1)*(dz(k) + dz(k - 1)) & + - q%sf(i, j, k - 1)*(dz(k) + dz(k + 1)) + dq%sf(i, j, k) = dq%sf(i, j, k)/ & + ((dz(k) + dz(k - 1))*(dz(k) + dz(k + 1))) end do end do - end if + end do end if end subroutine s_gradient_dir @@ -1465,7 +1923,7 @@ contains character(LEN=path_len + 2*name_len) :: file_loc logical :: file_exist - integer :: bub_id, tot_part, tot_part_wrtn, npart_wrtn + integer :: bub_id, tot_part integer :: i, k #ifdef MFC_MPI @@ -1476,6 +1934,9 @@ contains integer :: view integer, dimension(2) :: gsizes, lsizes, start_idx_part integer, dimension(num_procs) :: part_order, part_ord_mpi + integer, dimension(num_procs) :: proc_bubble_counts + real(wp), dimension(1:1, 1:lag_io_vars) :: dummy + dummy = 0._wp bub_id = 0._wp if (nBubs /= 0) then @@ -1488,78 +1949,60 @@ contains if (.not. parallel_io) return + lsizes(1) = bub_id + lsizes(2) = lag_io_vars + ! Total number of particles call MPI_ALLREDUCE(bub_id, tot_part, 1, MPI_integer, & MPI_SUM, MPI_COMM_WORLD, ierr) - ! Total number of particles written so far - call MPI_ALLREDUCE(npart_wrtn, tot_part_wrtn, 1, MPI_integer, & - MPI_SUM, MPI_COMM_WORLD, ierr) - - lsizes(1) = max(1, bub_id) - lsizes(2) = 21 - - ! if the particle number is zero, put 1 since MPI cannot deal with writing - ! zero particle - part_order(:) = 1 - part_order(proc_rank + 1) = max(1, bub_id) - - call MPI_ALLREDUCE(part_order, part_ord_mpi, num_procs, MPI_integer, & - MPI_MAX, MPI_COMM_WORLD, ierr) + call MPI_ALLGATHER(bub_id, 1, MPI_INTEGER, proc_bubble_counts, 1, MPI_INTEGER, & + MPI_COMM_WORLD, ierr) - gsizes(1) = sum(part_ord_mpi(1:num_procs)) - gsizes(2) = 21 - - start_idx_part(1) = sum(part_ord_mpi(1:proc_rank + 1)) - part_ord_mpi(proc_rank + 1) + ! Calculate starting index for this processor's particles + call MPI_EXSCAN(lsizes(1), start_idx_part(1), 1, MPI_INTEGER, MPI_SUM, MPI_COMM_WORLD, ierr) + if (proc_rank == 0) start_idx_part(1) = 0 start_idx_part(2) = 0 - write (file_loc, '(A,I0,A)') 'lag_bubbles_mpi_io_', t_step, '.dat' + gsizes(1) = tot_part + gsizes(2) = lag_io_vars + + write (file_loc, '(A,I0,A)') 'lag_bubbles_', t_step, '.dat' file_loc = trim(case_dir)//'/restart_data'//trim(mpiiofs)//trim(file_loc) - inquire (FILE=trim(file_loc), EXIST=file_exist) - if (file_exist .and. proc_rank == 0) then - call MPI_FILE_DELETE(file_loc, mpi_info_int, ierr) - end if - ! Writing down the total number of particles + ! Clean up existing file if (proc_rank == 0) then - open (9, FILE=trim(file_loc), FORM='unformatted', STATUS='unknown') - write (9) gsizes(1), mytime, dt - close (9) + inquire (FILE=trim(file_loc), EXIST=file_exist) + if (file_exist) then + call MPI_FILE_DELETE(file_loc, mpi_info_int, ierr) + end if end if - call MPI_type_CREATE_SUBARRAY(2, gsizes, lsizes, start_idx_part, & - MPI_ORDER_FORTRAN, mpi_p, view, ierr) - call MPI_type_COMMIT(view, ierr) + call MPI_BARRIER(MPI_COMM_WORLD, ierr) - allocate (MPI_IO_DATA_lag_bubbles(1:max(1, bub_id), 1:21)) - - ! Open the file to write all flow variables - write (file_loc, '(A,I0,A)') 'lag_bubbles_', t_step, '.dat' - file_loc = trim(case_dir)//'/restart_data'//trim(mpiiofs)//trim(file_loc) - inquire (FILE=trim(file_loc), EXIST=file_exist) - if (file_exist .and. proc_rank == 0) then - call MPI_FILE_DELETE(file_loc, mpi_info_int, ierr) - end if + if (proc_rank == 0) then + call MPI_FILE_OPEN(MPI_COMM_SELF, file_loc, & + ior(MPI_MODE_WRONLY, MPI_MODE_CREATE), & + mpi_info_int, ifile, ierr) - call MPI_FILE_OPEN(MPI_COMM_WORLD, file_loc, ior(MPI_MODE_WRONLY, MPI_MODE_CREATE), & - mpi_info_int, ifile, ierr) + ! Write header using MPI I/O for consistency + call MPI_FILE_WRITE(ifile, tot_part, 1, MPI_INTEGER, status, ierr) + call MPI_FILE_WRITE(ifile, mytime, 1, mpi_p, status, ierr) + call MPI_FILE_WRITE(ifile, dt, 1, mpi_p, status, ierr) + call MPI_FILE_WRITE(ifile, num_procs, 1, MPI_INTEGER, status, ierr) + call MPI_FILE_WRITE(ifile, proc_bubble_counts, num_procs, MPI_INTEGER, status, ierr) - disp = 0._wp + call MPI_FILE_CLOSE(ifile, ierr) + end if - call MPI_FILE_SET_VIEW(ifile, disp, mpi_p, view, & - 'native', mpi_info_null, ierr) + call MPI_BARRIER(MPI_COMM_WORLD, ierr) - ! Cycle through list - i = 1 - - if (bub_id == 0) then - MPI_IO_DATA_lag_bubbles(1, 1:21) = 0._wp - else + if (bub_id > 0) then + allocate (MPI_IO_DATA_lag_bubbles(max(1, bub_id), 1:lag_io_vars)) + i = 1 do k = 1, nBubs - if (particle_in_domain_physical(mtn_pos(k, 1:3, 1))) then - MPI_IO_DATA_lag_bubbles(i, 1) = real(lag_id(k, 1)) MPI_IO_DATA_lag_bubbles(i, 2:4) = mtn_pos(k, 1:3, 1) MPI_IO_DATA_lag_bubbles(i, 5:7) = mtn_posPrev(k, 1:3, 1) @@ -1575,21 +2018,47 @@ contains MPI_IO_DATA_lag_bubbles(i, 19) = gas_mg(k) MPI_IO_DATA_lag_bubbles(i, 20) = gas_betaT(k) MPI_IO_DATA_lag_bubbles(i, 21) = gas_betaC(k) - i = i + 1 - end if - end do - end if + call MPI_TYPE_CREATE_SUBARRAY(2, gsizes, lsizes, start_idx_part, & + MPI_ORDER_FORTRAN, mpi_p, view, ierr) + call MPI_TYPE_COMMIT(view, ierr) + + call MPI_FILE_OPEN(MPI_COMM_WORLD, file_loc, & + ior(MPI_MODE_WRONLY, MPI_MODE_CREATE), & + mpi_info_int, ifile, ierr) + + ! Skip header (written by rank 0) + disp = int(sizeof(tot_part) + 2*sizeof(mytime) + sizeof(num_procs) + & + num_procs*sizeof(proc_bubble_counts(1)), MPI_OFFSET_KIND) + call MPI_FILE_SET_VIEW(ifile, disp, mpi_p, view, 'native', mpi_info_int, ierr) + + call MPI_FILE_WRITE_ALL(ifile, MPI_IO_DATA_lag_bubbles, & + lag_io_vars*bub_id, mpi_p, status, ierr) - call MPI_FILE_write_ALL(ifile, MPI_IO_DATA_lag_bubbles, 21*max(1, bub_id), & - mpi_p, status, ierr) + call MPI_FILE_CLOSE(ifile, ierr) - call MPI_FILE_CLOSE(ifile, ierr) + deallocate (MPI_IO_DATA_lag_bubbles) + + else + call MPI_TYPE_CONTIGUOUS(0, mpi_p, view, ierr) + call MPI_TYPE_COMMIT(view, ierr) + + call MPI_FILE_OPEN(MPI_COMM_WORLD, file_loc, & + ior(MPI_MODE_WRONLY, MPI_MODE_CREATE), & + mpi_info_int, ifile, ierr) - deallocate (MPI_IO_DATA_lag_bubbles) + ! Skip header (written by rank 0) + disp = int(sizeof(tot_part) + 2*sizeof(mytime) + sizeof(num_procs) + & + num_procs*sizeof(proc_bubble_counts(1)), MPI_OFFSET_KIND) + call MPI_FILE_SET_VIEW(ifile, disp, mpi_p, view, 'native', mpi_info_int, ierr) + + call MPI_FILE_WRITE_ALL(ifile, dummy, 0, mpi_p, status, ierr) + + call MPI_FILE_CLOSE(ifile, ierr) + end if #endif @@ -1642,40 +2111,34 @@ contains !> The purpose of this subroutine is to remove one specific particle if dt is too small. !! @param bub_id Particle id - impure subroutine s_remove_lag_bubble(bub_id) - - integer, intent(in) :: bub_id - - integer :: i - - $:GPU_LOOP(parallelism='[seq]') - do i = bub_id, nBubs - 1 - lag_id(i, 1) = lag_id(i + 1, 1) - bub_R0(i) = bub_R0(i + 1) - Rmax_stats(i) = Rmax_stats(i + 1) - Rmin_stats(i) = Rmin_stats(i + 1) - gas_mg(i) = gas_mg(i + 1) - gas_betaT(i) = gas_betaT(i + 1) - gas_betaC(i) = gas_betaC(i + 1) - bub_dphidt(i) = bub_dphidt(i + 1) - gas_p(i, 1:2) = gas_p(i + 1, 1:2) - gas_mv(i, 1:2) = gas_mv(i + 1, 1:2) - intfc_rad(i, 1:2) = intfc_rad(i + 1, 1:2) - intfc_vel(i, 1:2) = intfc_vel(i + 1, 1:2) - mtn_pos(i, 1:3, 1:2) = mtn_pos(i + 1, 1:3, 1:2) - mtn_posPrev(i, 1:3, 1:2) = mtn_posPrev(i + 1, 1:3, 1:2) - mtn_vel(i, 1:3, 1:2) = mtn_vel(i + 1, 1:3, 1:2) - mtn_s(i, 1:3, 1:2) = mtn_s(i + 1, 1:3, 1:2) - intfc_draddt(i, 1:lag_num_ts) = intfc_draddt(i + 1, 1:lag_num_ts) - intfc_dveldt(i, 1:lag_num_ts) = intfc_dveldt(i + 1, 1:lag_num_ts) - gas_dpdt(i, 1:lag_num_ts) = gas_dpdt(i + 1, 1:lag_num_ts) - gas_dmvdt(i, 1:lag_num_ts) = gas_dmvdt(i + 1, 1:lag_num_ts) - end do - - nBubs = nBubs - 1 - $:GPU_UPDATE(device='[nBubs]') - - end subroutine s_remove_lag_bubble + impure subroutine s_copy_lag_bubble(dest, src) + + integer, intent(in) :: src, dest + + bub_R0(dest) = bub_R0(src) + Rmax_stats(dest) = Rmax_stats(src) + Rmin_stats(dest) = Rmin_stats(src) + gas_mg(dest) = gas_mg(src) + gas_betaT(dest) = gas_betaT(src) + gas_betaC(dest) = gas_betaC(src) + bub_dphidt(dest) = bub_dphidt(src) + lag_id(dest, 1) = lag_id(src, 1) + gas_p(dest, 1:2) = gas_p(src, 1:2) + gas_mv(dest, 1:2) = gas_mv(src, 1:2) + intfc_rad(dest, 1:2) = intfc_rad(src, 1:2) + intfc_vel(dest, 1:2) = intfc_vel(src, 1:2) + mtn_vel(dest, 1:3, 1:2) = mtn_vel(src, 1:3, 1:2) + mtn_s(dest, 1:3, 1:2) = mtn_s(src, 1:3, 1:2) + mtn_pos(dest, 1:3, 1:2) = mtn_pos(src, 1:3, 1:2) + mtn_posPrev(dest, 1:3, 1:2) = mtn_posPrev(src, 1:3, 1:2) + intfc_draddt(dest, 1:lag_num_ts) = intfc_draddt(src, 1:lag_num_ts) + intfc_dveldt(dest, 1:lag_num_ts) = intfc_dveldt(src, 1:lag_num_ts) + gas_dpdt(dest, 1:lag_num_ts) = gas_dpdt(src, 1:lag_num_ts) + gas_dmvdt(dest, 1:lag_num_ts) = gas_dmvdt(src, 1:lag_num_ts) + mtn_dposdt(dest, 1:3, 1:lag_num_ts) = mtn_dposdt(src, 1:3, 1:lag_num_ts) + mtn_dveldt(dest, 1:3, 1:lag_num_ts) = mtn_dveldt(src, 1:3, 1:lag_num_ts) + + end subroutine s_copy_lag_bubble !> The purpose of this subroutine is to deallocate variables impure subroutine s_finalize_lagrangian_solver() diff --git a/src/simulation/m_bubbles_EL_kernels.fpp b/src/simulation/m_bubbles_EL_kernels.fpp index 48ea3bad9a..484a433da0 100644 --- a/src/simulation/m_bubbles_EL_kernels.fpp +++ b/src/simulation/m_bubbles_EL_kernels.fpp @@ -128,7 +128,9 @@ contains s_coord(1:3) = lbk_s(l, 1:3, 2) center(1:2) = lbk_pos(l, 1:2, 2) if (p > 0) center(3) = lbk_pos(l, 3, 2) + cell = fd_number - buff_size call s_get_cell(s_coord, cell) + !print*, s_coord call s_compute_stddsv(cell, volpart, stddsv) strength_vol = volpart @@ -248,7 +250,7 @@ contains theta = 0._wp Nr = ceiling(lag_params%charwidth/(y_cb(cellaux(2)) - y_cb(cellaux(2) - 1))) Nr_count = 1._wp - mapCells*1._wp - dzp = y_cb(cellaux(2) + 1) - y_cb(cellaux(2)) + dzp = y_cb(cellaux(2)) - y_cb(cellaux(2) - 1) Lz2 = (center(3) - (dzp*(0.5_wp + Nr_count) - lag_params%charwidth/2._wp))**2._wp distance = sqrt((center(1) - nodecoord(1))**2._wp + (center(2) - nodecoord(2))**2._wp + Lz2) func = dzp/lag_params%charwidth*exp(-0.5_wp*(distance/stddsv)**2._wp)/(sqrt(2._wp*pi)*stddsv)**3._wp @@ -277,21 +279,27 @@ contains celloutside = .false. if (num_dims == 2) then - if ((cellaux(1) < -buff_size) .or. (cellaux(2) < -buff_size)) then + if ((cellaux(1) < fd_number - buff_size) .or. & + (cellaux(2) < fd_number - buff_size)) then celloutside = .true. end if - if (cyl_coord .and. y_cc(cellaux(2)) < 0._wp) then + if (cyl_coord .and. cellaux(2) < 0) then celloutside = .true. end if - if ((cellaux(2) > n + buff_size) .or. (cellaux(1) > m + buff_size)) then + if ((cellaux(2) > n + buff_size - fd_number) .or. & + (cellaux(1) > m + buff_size - fd_number)) then celloutside = .true. end if else - if ((cellaux(3) < -buff_size) .or. (cellaux(1) < -buff_size) .or. (cellaux(2) < -buff_size)) then + if ((cellaux(3) < fd_number - buff_size) .or. & + (cellaux(1) < fd_number - buff_size) .or. & + (cellaux(2) < fd_number - buff_size)) then celloutside = .true. end if - if ((cellaux(3) > p + buff_size) .or. (cellaux(2) > n + buff_size) .or. (cellaux(1) > m + buff_size)) then + if ((cellaux(3) > p + buff_size - fd_number) .or. & + (cellaux(2) > n + buff_size - fd_number) .or. & + (cellaux(1) > m + buff_size - fd_number)) then celloutside = .true. end if end if @@ -417,4 +425,231 @@ contains end subroutine s_get_cell + !! This function interpolates the velocity of Eulerian field at the position + !! of the bubble. + !! @param pos Position of the bubble in directiion i + !! @param cell Computational coordinates of the bubble + !! @param i Direction of the velocity (1: x, 2: y, 3: z) + !! @param q_prim_vf Eulerian field with primitive variables + !! @return v Interpolated velocity at the position of the bubble + pure function f_interpolate_velocity(pos, cell, i, q_prim_vf) result(v) + $:GPU_ROUTINE(parallelism='[seq]') + real(wp), intent(in) :: pos + integer, dimension(3), intent(in) :: cell + integer, intent(in) :: i + type(scalar_field), dimension(sys_size), intent(in) :: q_prim_vf + + real(wp) :: v + real(wp), dimension(fd_order + 1) :: xi, eta, L + + if (fd_order == 2) then + if (i == 1) then + xi(1) = x_cc(cell(1) - 1) + eta(1) = q_prim_vf(momxb)%sf(cell(1) - 1, cell(2), cell(3)) + xi(2) = x_cc(cell(1)) + eta(2) = q_prim_vf(momxb)%sf(cell(1), cell(2), cell(3)) + xi(3) = x_cc(cell(1) + 1) + eta(3) = q_prim_vf(momxb)%sf(cell(1) + 1, cell(2), cell(3)) + elseif (i == 2) then + xi(1) = y_cc(cell(2) - 1) + eta(1) = q_prim_vf(momxb + 1)%sf(cell(1), cell(2) - 1, cell(3)) + xi(2) = y_cc(cell(2)) + eta(2) = q_prim_vf(momxb + 1)%sf(cell(1), cell(2), cell(3)) + xi(3) = y_cc(cell(2) + 1) + eta(3) = q_prim_vf(momxb + 1)%sf(cell(1), cell(2) + 1, cell(3)) + elseif (i == 3) then + xi(1) = z_cc(cell(3) - 1) + eta(1) = q_prim_vf(momxe)%sf(cell(1), cell(2), cell(3) - 1) + xi(2) = z_cc(cell(3)) + eta(2) = q_prim_vf(momxe)%sf(cell(1), cell(2), cell(3)) + xi(3) = z_cc(cell(3) + 1) + eta(3) = q_prim_vf(momxe)%sf(cell(1), cell(2), cell(3) + 1) + end if + + L(1) = ((pos - xi(2))*(pos - xi(3)))/((xi(1) - xi(2))*(xi(1) - xi(3))) + L(2) = ((pos - xi(1))*(pos - xi(3)))/((xi(2) - xi(1))*(xi(2) - xi(3))) + L(3) = ((pos - xi(1))*(pos - xi(2)))/((xi(3) - xi(1))*(xi(3) - xi(2))) + + v = L(1)*eta(1) + L(2)*eta(2) + L(3)*eta(3) + elseif (fd_order == 4) then + if (i == 1) then + xi(1) = x_cc(cell(1) - 2) + eta(1) = q_prim_vf(momxb)%sf(cell(1) - 2, cell(2), cell(3)) + xi(2) = x_cc(cell(1) - 1) + eta(2) = q_prim_vf(momxb)%sf(cell(1) - 1, cell(2), cell(3)) + xi(3) = x_cc(cell(1)) + eta(3) = q_prim_vf(momxb)%sf(cell(1), cell(2), cell(3)) + xi(4) = x_cc(cell(1) + 1) + eta(4) = q_prim_vf(momxb)%sf(cell(1) + 1, cell(2), cell(3)) + xi(5) = x_cc(cell(1) + 2) + eta(5) = q_prim_vf(momxb)%sf(cell(1) + 2, cell(2), cell(3)) + elseif (i == 2) then + xi(1) = y_cc(cell(2) - 2) + eta(1) = q_prim_vf(momxb + 1)%sf(cell(1), cell(2) - 2, cell(3)) + xi(2) = y_cc(cell(2) - 1) + eta(2) = q_prim_vf(momxb + 1)%sf(cell(1), cell(2) - 1, cell(3)) + xi(3) = y_cc(cell(2)) + eta(3) = q_prim_vf(momxb + 1)%sf(cell(1), cell(2), cell(3)) + xi(4) = y_cc(cell(2) + 1) + eta(4) = q_prim_vf(momxb + 1)%sf(cell(1), cell(2) + 1, cell(3)) + xi(5) = y_cc(cell(2) + 2) + eta(5) = q_prim_vf(momxb + 1)%sf(cell(1), cell(2) + 2, cell(3)) + elseif (i == 3) then + xi(1) = z_cc(cell(3) - 2) + eta(1) = q_prim_vf(momxe)%sf(cell(1), cell(2), cell(3) - 2) + xi(2) = z_cc(cell(3) - 1) + eta(2) = q_prim_vf(momxe)%sf(cell(1), cell(2), cell(3) - 1) + xi(3) = z_cc(cell(3)) + eta(3) = q_prim_vf(momxe)%sf(cell(1), cell(2), cell(3)) + xi(4) = z_cc(cell(3) + 1) + eta(4) = q_prim_vf(momxe)%sf(cell(1), cell(2), cell(3) + 1) + xi(5) = z_cc(cell(3) + 2) + eta(5) = q_prim_vf(momxe)%sf(cell(1), cell(2), cell(3) + 2) + end if + + L(1) = ((pos - xi(2))*(pos - xi(3))*(pos - xi(4))*(pos - xi(5)))/ & + ((xi(1) - xi(2))*(xi(1) - xi(3))*(xi(1) - xi(3))*(xi(2) - xi(5))) + L(2) = ((pos - xi(1))*(pos - xi(3))*(pos - xi(4))*(pos - xi(5)))/ & + ((xi(2) - xi(1))*(xi(2) - xi(3))*(xi(2) - xi(3))*(xi(2) - xi(5))) + L(3) = ((pos - xi(1))*(pos - xi(2))*(pos - xi(4))*(pos - xi(5)))/ & + ((xi(3) - xi(1))*(xi(3) - xi(2))*(xi(3) - xi(4))*(xi(3) - xi(5))) + L(4) = ((pos - xi(1))*(pos - xi(2))*(pos - xi(3))*(pos - xi(4)))/ & + ((xi(4) - xi(1))*(xi(4) - xi(2))*(xi(4) - xi(3))*(xi(4) - xi(5))) + L(5) = ((pos - xi(1))*(pos - xi(2))*(pos - xi(3))*(pos - xi(4)))/ & + ((xi(5) - xi(1))*(xi(5) - xi(2))*(xi(5) - xi(3))*(xi(5) - xi(4))) + + v = L(1)*eta(1) + L(2)*eta(2) + L(3)*eta(3) + L(4)*eta(4) + L(5)*eta(5) + end if + + end function f_interpolate_velocity + + !! This function calculates the force on a bubble + !! based on the pressure gradient, velocity, and drag model. + !! @param pos Position of the bubble in direction i + !! @param rad Radius of the bubble + !! @param rdot Radial velocity of the bubble + !! @param vel Velocity of the bubble + !! @param mg Mass of the gas in the bubble + !! @param mv Mass of the liquid in the bubble + !! @param Re Reynolds number + !! @param rho Density of the fluid + !! @param cell Computational coordinates of the bubble + !! @param i Direction of the velocity (1: x, 2: y, 3: z) + !! @param q_prim_vf Eulerian field with primitive variables + !! @return a Acceleration of the bubble in direction i + pure function f_get_bubble_force(pos, rad, rdot, vel, mg, mv, Re, rho, cell, i, q_prim_vf) result(force) + $:GPU_ROUTINE(parallelism='[seq]') + real(wp), intent(in) :: pos, rad, rdot, mg, mv, Re, rho + real(wp), intent(in), dimension(3) :: vel + integer, dimension(3), intent(in) :: cell + integer, intent(in) :: i + type(scalar_field), dimension(sys_size), intent(in) :: q_prim_vf + + real(wp) :: a, dp, vol, force, Re_b, C_d, v_rel_mag, rho_b + real(wp), dimension(3) :: v_rel + real(wp), dimension(fd_order - 1) :: xi, eta, L + integer :: j + + if (fd_order == 2) then + if (i == 1) then + dp = (q_prim_vf(E_idx)%sf(cell(1) + 1, cell(2), cell(3)) - & + q_prim_vf(E_idx)%sf(cell(1) - 1, cell(2), cell(3)))/ & + (x_cc(cell(1) + 1) - x_cc(cell(1) - 1)) + elseif (i == 2) then + dp = (q_prim_vf(E_idx)%sf(cell(1), cell(2) + 1, cell(3)) - & + q_prim_vf(E_idx)%sf(cell(1), cell(2) - 1, cell(3)))/ & + (y_cc(cell(2) + 1) - y_cc(cell(2) - 1)) + elseif (i == 3) then + dp = (q_prim_vf(E_idx)%sf(cell(1), cell(2), cell(3) + 1) - & + q_prim_vf(E_idx)%sf(cell(1), cell(2), cell(3) - 1))/ & + (z_cc(cell(3) + 1) - z_cc(cell(3) - 1)) + end if + elseif (fd_order == 4) then + if (i == 1) then + xi(1) = x_cc(cell(1) - 1) + eta(1) = (q_prim_vf(E_idx)%sf(cell(1), cell(2), cell(3)) - & + q_prim_vf(E_idx)%sf(cell(1) - 2, cell(2), cell(3)))/ & + (x_cc(cell(1)) - x_cc(cell(1) - 2)) + xi(2) = x_cc(cell(1)) + eta(2) = (q_prim_vf(E_idx)%sf(cell(1) + 1, cell(2), cell(3)) - & + q_prim_vf(E_idx)%sf(cell(1) - 1, cell(2), cell(3)))/ & + (x_cc(cell(1) + 1) - x_cc(cell(1) - 1)) + xi(3) = x_cc(cell(1) + 1) + eta(3) = (q_prim_vf(E_idx)%sf(cell(1) + 2, cell(2), cell(3)) - & + q_prim_vf(E_idx)%sf(cell(1), cell(2), cell(3)))/ & + (x_cc(cell(1) + 2) - x_cc(cell(1))) + elseif (i == 2) then + xi(1) = y_cc(cell(2) - 1) + eta(1) = (q_prim_vf(E_idx)%sf(cell(1), cell(2), cell(3)) - & + q_prim_vf(E_idx)%sf(cell(1), cell(2) - 2, cell(3)))/ & + (y_cc(cell(2)) - y_cc(cell(2) - 2)) + xi(2) = y_cc(cell(2)) + eta(2) = (q_prim_vf(E_idx)%sf(cell(1), cell(2) + 1, cell(3)) - & + q_prim_vf(E_idx)%sf(cell(1), cell(2) - 1, cell(3)))/ & + (y_cc(cell(2) + 1) - y_cc(cell(2) - 1)) + xi(3) = y_cc(cell(2) + 1) + eta(3) = (q_prim_vf(E_idx)%sf(cell(1), cell(2) + 2, cell(3)) - & + q_prim_vf(E_idx)%sf(cell(1), cell(2), cell(3)))/ & + (y_cc(cell(2) + 2) - y_cc(cell(2))) + elseif (i == 3) then + xi(1) = z_cc(cell(3) - 1) + eta(1) = (q_prim_vf(E_idx)%sf(cell(1), cell(2), cell(3)) - & + q_prim_vf(E_idx)%sf(cell(1), cell(2), cell(3) - 2))/ & + (z_cc(cell(3)) - z_cc(cell(3) - 2)) + xi(2) = z_cc(cell(3)) + eta(2) = (q_prim_vf(E_idx)%sf(cell(1), cell(2), cell(3) + 1) - & + q_prim_vf(E_idx)%sf(cell(1), cell(2), cell(3) - 1))/ & + (z_cc(cell(3) + 1) - z_cc(cell(3) - 1)) + xi(3) = z_cc(cell(3) + 1) + eta(3) = (q_prim_vf(E_idx)%sf(cell(1), cell(2), cell(3) + 2) - & + q_prim_vf(E_idx)%sf(cell(1), cell(2), cell(3)))/ & + (z_cc(cell(3) + 2) - z_cc(cell(3))) + end if + + L(1) = ((pos - xi(2))*(pos - xi(3)))/((xi(1) - xi(2))*(xi(1) - xi(3))) + L(2) = ((pos - xi(1))*(pos - xi(3)))/((xi(2) - xi(1))*(xi(2) - xi(3))) + L(3) = ((pos - xi(1))*(pos - xi(2)))/((xi(3) - xi(1))*(xi(3) - xi(2))) + + dp = L(1)*eta(1) + L(2)*eta(2) + L(3)*eta(3) + end if + + vol = (4._wp/3._wp) * pi * (rad**3._wp) + v_rel_mag = 0._wp + + do j = 1, num_dims + v_rel(j) = vel(j) - f_interpolate_velocity(pos, cell, j, q_prim_vf) + v_rel_mag = v_rel_mag + v_rel(j)**2._wp + end do + + force = 0._wp + + if (lag_params%drag_model == 1) then ! Free slip Stokes drag + force = force - (4._wp*pi*rad*v_rel(i))/Re + else if (lag_params%drag_model == 2) then ! No slip Stokes drag + force = force - (6._wp*pi*rad*v_rel(i))/Re + else if (lag_params%drag_model == 3) then ! Levich drag + force = force - (12._wp*pi*rad*v_rel(i))/Re + else if (lag_params%drag_model > 0) then ! Drag coefficient model + v_rel_mag = sqrt(v_rel_mag) + Re_b = max(1e-3, rho * v_rel_mag * 2._wp * rad * Re) + if (lag_params%drag_model == 4) then ! Mei et al. 1994 + C_d = 16._wp / Re_b * (1._wp + (8._wp / Re_b + 0.5_wp * (1._wp + 3.315_wp * Re_b**(-0.5_wp))) ** -1._wp) + end if + force = force - 0.5_wp * C_d * rho * pi * rad**2._wp * v_rel(i) * v_rel_mag + end if + + if (lag_params%pressure_force) then + force = force - vol * dp + end if + + if (lag_params%gravity_force) then + force = force + (mg + mv) * accel_bf(i) + end if + + if (lag_params%momentum_transfer_force) then + force = force - 4._wp * pi * rho * rad**2._wp * v_rel(i) * rdot + end if + + end function f_get_bubble_force + end module m_bubbles_EL_kernels diff --git a/src/simulation/m_data_output.fpp b/src/simulation/m_data_output.fpp index 5f5b79c481..40b80df48c 100644 --- a/src/simulation/m_data_output.fpp +++ b/src/simulation/m_data_output.fpp @@ -36,7 +36,7 @@ module m_data_output implicit none - private; + private; public :: s_initialize_data_output_module, & s_open_run_time_information_file, & s_open_com_files, & @@ -157,24 +157,16 @@ contains write (3, '(A)') ''; write (3, '(A)') '' ! Generating table header for the stability criteria to be outputted - if (cfl_dt) then - if (viscous) then - write (3, '(A)') ' Time-steps dt = Time ICFL '// & - 'Max VCFL Max Rc Min =' - else - write (3, '(A)') ' Time-steps dt Time '// & - ' ICFL Max ' - end if - else - if (viscous) then - write (3, '(A)') ' Time-steps Time ICFL '// & - 'Max VCFL Max Rc Min ' - else - write (3, '(A)') ' Time-steps Time '// & - ' ICFL Max ' - end if + write (3, '(13X,A9,13X,A10,13X,A10,13X,A10)', advance="no") & + trim('Time-step'), trim('dt'), trim('Time'), trim('ICFL Max') + + if (viscous) then + write (3, '(13X,A10,13X,A16)', advance="no") & + trim('VCFL Max'), trim('Rc Min') end if + write (3, *) ! new line + end subroutine s_open_run_time_information_file !> This opens a formatted data file where the root processor @@ -296,7 +288,7 @@ contains call s_compute_speed_of_sound(pres, rho, gamma, pi_inf, H, alpha, vel_sum, 0._wp, c) if (viscous) then - call s_compute_stability_from_dt(vel, c, rho, Re, j, k, l, icfl_sf, vcfl_sf, Rc_sf) + call s_compute_stability_from_dt(vel, c, rho, Re, j, k, l, icfl_sf, vcfl_sf=vcfl_sf, Rc_sf=Rc_sf) else call s_compute_stability_from_dt(vel, c, rho, Re, j, k, l, icfl_sf) end if @@ -311,14 +303,10 @@ contains #ifdef _CRAYFTN $:GPU_UPDATE(host='[icfl_sf]') - - if (viscous) then - $:GPU_UPDATE(host='[vcfl_sf,Rc_sf]') - end if - icfl_max_loc = maxval(icfl_sf) if (viscous) then + $:GPU_UPDATE(host='[vcfl_sf,Rc_sf]') vcfl_max_loc = maxval(vcfl_sf) Rc_min_loc = minval(Rc_sf) end if @@ -358,16 +346,17 @@ contains ! Outputting global stability criteria extrema at current time-step if (proc_rank == 0) then + write (3, '(13X,I9,13X,F10.6,13X,F10.6,13X,F10.6)', advance="no") & + t_step, dt, mytime, icfl_max_glb + if (viscous) then - write (3, '(6X,I8,F10.6,6X,6X,F10.6,6X,F9.6,6X,F9.6,6X,F10.6)') & - t_step, dt, t_step*dt, icfl_max_glb, & + write (3, '(13X,F10.6,13X,ES16.6)', advance="no") & vcfl_max_glb, & Rc_min_glb - else - write (3, '(13X,I8,14X,F10.6,14X,F10.6,13X,F9.6)') & - t_step, dt, t_step*dt, icfl_max_glb end if + write (3, *) ! new line + if (.not. f_approx_equal(icfl_max_glb, icfl_max_glb)) then call s_mpi_abort('ICFL is NaN. Exiting.') elseif (icfl_max_glb > 1._wp) then @@ -383,6 +372,7 @@ contains call s_mpi_abort('VCFL is greater than 1.0. Exiting.') end if end if + end if call s_mpi_barrier() @@ -1756,7 +1746,7 @@ contains write (3, '(A,F9.6)') 'ICFL Max: ', icfl_max if (viscous) write (3, '(A,F9.6)') 'VCFL Max: ', vcfl_max - if (viscous) write (3, '(A,F10.6)') 'Rc Min: ', Rc_min + if (viscous) write (3, '(A,ES16.6)') 'Rc Min: ', Rc_min call cpu_time(run_time) @@ -1805,7 +1795,7 @@ contains @:ALLOCATE(Rc_sf (0:m, 0:n, 0:p)) vcfl_max = 0._wp - Rc_min = 1.e3_wp + Rc_min = 1.e12_wp end if end if diff --git a/src/simulation/m_global_parameters.fpp b/src/simulation/m_global_parameters.fpp index 37075a3c4d..a101ecb380 100644 --- a/src/simulation/m_global_parameters.fpp +++ b/src/simulation/m_global_parameters.fpp @@ -63,6 +63,7 @@ module m_global_parameters !> @{ real(wp), target, allocatable, dimension(:) :: x_cb, y_cb, z_cb + type(bounds_info), dimension(3) :: glb_bounds !< !> @} !> @name Cell-center (CC) locations in the x-, y- and z-directions, respectively @@ -70,7 +71,6 @@ module m_global_parameters real(wp), target, allocatable, dimension(:) :: x_cc, y_cc, z_cc !> @} - !type(bounds_info) :: x_domain, y_domain, z_domain !< !! Locations of the domain bounds in the x-, y- and z-coordinate directions !> @name Cell-width distributions in the x-, y- and z-directions, respectively !> @{ @@ -231,6 +231,7 @@ module m_global_parameters integer :: num_bc_patches logical :: bc_io + logical, dimension(3) :: periodic_bc !> @name Boundary conditions (BC) in the x-, y- and z-directions, respectively !> @{ type(int_bounds_info) :: bc_x, bc_y, bc_z @@ -239,10 +240,6 @@ module m_global_parameters $:GPU_DECLARE(create='[bc_y%vb1, bc_y%vb2, bc_y%vb3, bc_y%ve1, bc_y%ve2, bc_y%ve3]') $:GPU_DECLARE(create='[bc_z%vb1, bc_z%vb2, bc_z%vb3, bc_z%ve1, bc_z%ve2, bc_z%ve3]') - type(bounds_info) :: x_domain, y_domain, z_domain - real(wp) :: x_a, y_a, z_a - real(wp) :: x_b, y_b, z_b - logical :: parallel_io !< Format of the data files logical :: file_per_process !< shared file or not when using parallel io integer :: precision !< Precision of output files @@ -252,6 +249,15 @@ module m_global_parameters integer, allocatable, dimension(:) :: proc_coords !< !! Processor coordinates in MPI_CART_COMM + type(bounds_info), allocatable, dimension(:) :: pcomm_coords + $:GPU_DECLARE(create='[pcomm_coords]') + !! Coordinates for EL particle transfer + + type(int_bounds_info), dimension(3) :: nidx !< Indices for neighboring processors + + integer, allocatable, dimension(:, :, :) :: neighbor_ranks + !! Neighbor ranks for lagrangian particle communication + integer, allocatable, dimension(:) :: start_idx !< !! Starting cell-center index of local processor in global grid @@ -647,6 +653,7 @@ contains num_bc_patches = 0 bc_io = .false. + periodic_bc = .false. bc_x%beg = dflt_int; bc_x%end = dflt_int bc_y%beg = dflt_int; bc_y%end = dflt_int @@ -659,9 +666,9 @@ contains #:endfor #:endfor - x_domain%beg = dflt_int; x_domain%end = dflt_int - y_domain%beg = dflt_int; y_domain%end = dflt_int - z_domain%beg = dflt_int; z_domain%end = dflt_int + glb_bounds(1)%beg = dflt_int; glb_bounds(1)%end = dflt_int + glb_bounds(2)%beg = dflt_int; glb_bounds(2)%end = dflt_int + glb_bounds(3)%beg = dflt_int; glb_bounds(3)%end = dflt_int ! Fluids physical parameters do i = 1, num_fluids_max @@ -807,6 +814,12 @@ contains lag_params%write_bubbles = .false. lag_params%write_bubbles_stats = .false. lag_params%nBubs_glb = dflt_int + lag_params%vel_model = dflt_int + lag_params%drag_model = dflt_int + lag_params%pressure_force = .true. + lag_params%gravity_force = .false. + lag_params%momentum_transfer_force = .false. + lag_params%c_d = dflt_real lag_params%epsilonb = 1._wp lag_params%charwidth = dflt_real lag_params%valmaxvoid = dflt_real @@ -1250,11 +1263,15 @@ contains fd_number = max(1, fd_order/2) end if + if (bubbles_lagrange) then + fd_number = max(1, fd_order/2) + end if + call s_configure_coordinate_bounds(recon_type, weno_polyn, muscl_polyn, & igr_order, buff_size, & idwint, idwbuff, viscous, & bubbles_lagrange, m, n, p, & - num_dims, igr, ib) + num_dims, igr, ib, fd_number) $:GPU_UPDATE(device='[idwint, idwbuff]') ! Configuring Coordinate Direction Indexes @@ -1378,6 +1395,7 @@ contains #:endif allocate (proc_coords(1:num_dims)) + @:ALLOCATE(pcomm_coords(1:num_dims)) if (parallel_io .neqv. .true.) return @@ -1414,6 +1432,7 @@ contains end if deallocate (proc_coords) + @:DEALLOCATE(pcomm_coords) if (parallel_io) then deallocate (start_idx) diff --git a/src/simulation/m_ibm.fpp b/src/simulation/m_ibm.fpp index 6fe4932574..0fd6869ba6 100644 --- a/src/simulation/m_ibm.fpp +++ b/src/simulation/m_ibm.fpp @@ -546,17 +546,17 @@ contains ghost_points_in(count)%slip = patch_ib(patch_id)%slip ! ghost_points(count)%rank = proc_rank - if ((x_cc(i) - dx(i)) < x_domain%beg) then + if ((x_cc(i) - dx(i)) < glb_bounds(1)%beg) then ghost_points_in(count)%DB(1) = -1 - else if ((x_cc(i) + dx(i)) > x_domain%end) then + else if ((x_cc(i) + dx(i)) > glb_bounds(1)%end) then ghost_points_in(count)%DB(1) = 1 else ghost_points_in(count)%DB(1) = 0 end if - if ((y_cc(j) - dy(j)) < y_domain%beg) then + if ((y_cc(j) - dy(j)) < glb_bounds(2)%beg) then ghost_points_in(count)%DB(2) = -1 - else if ((y_cc(j) + dy(j)) > y_domain%end) then + else if ((y_cc(j) + dy(j)) > glb_bounds(2)%end) then ghost_points_in(count)%DB(2) = 1 else ghost_points_in(count)%DB(2) = 0 @@ -589,25 +589,25 @@ contains ib_markers%sf(i, j, k) ghost_points_in(count)%slip = patch_ib(patch_id)%slip - if ((x_cc(i) - dx(i)) < x_domain%beg) then + if ((x_cc(i) - dx(i)) < glb_bounds(1)%beg) then ghost_points_in(count)%DB(1) = -1 - else if ((x_cc(i) + dx(i)) > x_domain%end) then + else if ((x_cc(i) + dx(i)) > glb_bounds(1)%end) then ghost_points_in(count)%DB(1) = 1 else ghost_points_in(count)%DB(1) = 0 end if - if ((y_cc(j) - dy(j)) < y_domain%beg) then + if ((y_cc(j) - dy(j)) < glb_bounds(2)%beg) then ghost_points_in(count)%DB(2) = -1 - else if ((y_cc(j) + dy(j)) > y_domain%end) then + else if ((y_cc(j) + dy(j)) > glb_bounds(2)%end) then ghost_points_in(count)%DB(2) = 1 else ghost_points_in(count)%DB(2) = 0 end if - if ((z_cc(k) - dz(k)) < z_domain%beg) then + if ((z_cc(k) - dz(k)) < glb_bounds(3)%beg) then ghost_points_in(count)%DB(3) = -1 - else if ((z_cc(k) + dz(k)) > z_domain%end) then + else if ((z_cc(k) + dz(k)) > glb_bounds(3)%end) then ghost_points_in(count)%DB(3) = 1 else ghost_points_in(count)%DB(3) = 0 diff --git a/src/simulation/m_mpi_proxy.fpp b/src/simulation/m_mpi_proxy.fpp index 81460b5c70..6c27b8cc7d 100644 --- a/src/simulation/m_mpi_proxy.fpp +++ b/src/simulation/m_mpi_proxy.fpp @@ -32,6 +32,8 @@ module m_mpi_proxy implicit none + integer, parameter :: n_neighbor = 26 + integer, private, allocatable, dimension(:) :: ib_buff_send !< !! This variable is utilized to pack and send the buffer of the immersed !! boundary markers, for a single computational domain boundary at the @@ -45,10 +47,26 @@ module m_mpi_proxy integer :: i_halo_size $:GPU_DECLARE(create='[i_halo_size]') + integer, dimension(-1:1, -1:1, -1:1) :: p_send_counts, p_recv_counts + integer, dimension(:, :, :, :), allocatable :: p_send_ids + character(len=1), dimension(:), allocatable :: p_send_buff, p_recv_buff + integer :: p_buff_size, p_var_size + + !! EL Bubbles communication variables + integer, parameter :: MAX_NEIGHBORS = 27 + integer :: send_requests(MAX_NEIGHBORS), recv_requests(MAX_NEIGHBORS) + integer :: recv_offsets(MAX_NEIGHBORS) + integer :: neighbor_list(MAX_NEIGHBORS, 3) + integer :: n_neighbors + + $:GPU_DECLARE(create='[p_send_counts]') + contains subroutine s_initialize_mpi_proxy_module() + integer :: i, j, k + #ifdef MFC_MPI if (ib) then if (n > 0) then @@ -73,6 +91,44 @@ contains end subroutine s_initialize_mpi_proxy_module + !! This subroutine initializes the MPI buffers and variables + !! required for the particle communication. + !! @param lag_num_ts Number of stages in time-stepping scheme + subroutine s_initialize_particles_mpi(lag_num_ts) + + integer :: i, j, k + integer :: real_size, int_size, nReal, lag_num_ts + integer :: ierr !< Generic flag used to identify and report MPI errors + +#ifdef MFC_MPI + call MPI_Pack_size(1, mpi_p, MPI_COMM_WORLD, real_size, ierr) + call MPI_Pack_size(1, MPI_INTEGER, MPI_COMM_WORLD, int_size, ierr) + + nReal = 7 + 16*2 + 10*lag_num_ts + p_var_size = (nReal*real_size + int_size) + p_buff_size = lag_params%nBubs_glb*p_var_size + + @:ALLOCATE(p_send_buff(0:p_buff_size), p_recv_buff(0:p_buff_size)) + @:ALLOCATE(p_send_ids(nidx(1)%beg:nidx(1)%end, nidx(2)%beg:nidx(2)%end, nidx(3)%beg:nidx(3)%end, 0:lag_params%nBubs_glb)) + + ! First, collect all neighbor information + n_neighbors = 0 + do k = nidx(3)%beg, nidx(3)%end + do j = nidx(2)%beg, nidx(2)%end + do i = nidx(1)%beg, nidx(1)%end + if (abs(i) + abs(j) + abs(k) /= 0) then + n_neighbors = n_neighbors + 1 + neighbor_list(n_neighbors, 1) = i + neighbor_list(n_neighbors, 2) = j + neighbor_list(n_neighbors, 3) = k + end if + end do + end do + end do +#endif + + end subroutine s_initialize_particles_mpi + !> Since only the processor with rank 0 reads and verifies !! the consistency of user inputs, these are initially not !! available to the other processors. Then, the purpose of @@ -100,7 +156,7 @@ contains & 'num_probes', 'num_integrals', 'bubble_model', 'thermal', & & 'num_source', 'relax_model', 'num_ibs', 'n_start', & & 'num_bc_patches', 'num_igr_iters', 'num_igr_warm_start_iters', & - & 'adap_dt_max_iters' ] + & 'adap_dt_max_iters', 'lag_params%vel_model' ] call MPI_BCAST(${VAR}$, 1, MPI_INTEGER, 0, MPI_COMM_WORLD, ierr) #:endfor @@ -132,7 +188,8 @@ contains if (bubbles_lagrange) then #:for VAR in [ 'heatTransfer_model', 'massTransfer_model', 'pressure_corrector', & - & 'write_bubbles', 'write_bubbles_stats'] + & 'write_bubbles', 'write_bubbles_stats', 'vel_model', 'drag_model', & + & 'pressure_force', 'gravity_force', 'momentum_transfer_force'] call MPI_BCAST(lag_params%${VAR}$, 1, MPI_LOGICAL, 0, MPI_COMM_WORLD, ierr) #:endfor @@ -152,11 +209,8 @@ contains & 'bc_y%vb1','bc_y%vb2','bc_y%vb3','bc_y%ve1','bc_y%ve2','bc_y%ve3', & & 'bc_z%vb1','bc_z%vb2','bc_z%vb3','bc_z%ve1','bc_z%ve2','bc_z%ve3', & & 'bc_x%pres_in','bc_x%pres_out','bc_y%pres_in','bc_y%pres_out', 'bc_z%pres_in','bc_z%pres_out', & - & 'x_domain%beg', 'x_domain%end', 'y_domain%beg', 'y_domain%end', & - & 'z_domain%beg', 'z_domain%end', 'x_a', 'x_b', 'y_a', 'y_b', 'z_a', & - & 'z_b', 't_stop', 't_save', 'cfl_target', 'Bx0', 'alf_factor', & - & 'tau_star', 'cont_damage_s', 'alpha_bar', 'adap_dt_tol', & - & 'ic_eps', 'ic_beta' ] + & 't_stop', 't_save', 'cfl_target', 'Bx0', 'alf_factor', 'tau_star', & + & 'cont_damage_s', 'alpha_bar', 'adap_dt_tol', 'ic_eps', 'ic_beta' ] call MPI_BCAST(${VAR}$, 1, mpi_p, 0, MPI_COMM_WORLD, ierr) #:endfor @@ -246,6 +300,9 @@ contains call MPI_BCAST(nv_uvm_igr_temps_on_gpu, 1, MPI_INTEGER, 0, MPI_COMM_WORLD, ierr) call MPI_BCAST(nv_uvm_pref_gpu, 1, MPI_LOGICAL, 0, MPI_COMM_WORLD, ierr) + ! Extra BC Variable + call MPI_BCAST(periodic_bc, 3, MPI_LOGICAL, 0, MPI_COMM_WORLD, ierr) + #endif end subroutine s_mpi_bcast_user_inputs @@ -428,6 +485,489 @@ contains end subroutine s_mpi_sendrecv_ib_buffers + !> This subroutine adds particles to the transfer list for the MPI + !! communication. + !! @param nBub Current LOCAL number of bubbles + !! @param pos Current position of each bubble + !! @param posPrev Previous position of each bubble (optional, not used + !! for communication of initial condition) + impure subroutine s_add_particles_to_transfer_list(nBub, pos, posPrev) + + real(wp), dimension(:, :) :: pos + real(wp), dimension(:, :), optional :: posPrev + integer :: bubID, nbub + integer :: i, j, k + + do k = nidx(3)%beg, nidx(3)%end + do j = nidx(2)%beg, nidx(2)%end + do i = nidx(1)%beg, nidx(1)%end + p_send_counts(i, j, k) = 0 + end do + end do + end do + + do k = 1, nbub + if (f_crosses_boundary(k, 1, -1, pos, posPrev)) then + call s_add_particle_to_direction(k, -1, 0, 0) + if (n > 0) then + if (f_crosses_boundary(k, 2, -1, pos, posPrev)) then + call s_add_particle_to_direction(k, -1, -1, 0) + call s_add_particle_to_direction(k, 0, -1, 0) + if (p > 0) then + if (f_crosses_boundary(k, 3, -1, pos, posPrev)) then + call s_add_particle_to_direction(k, -1, -1, -1) + call s_add_particle_to_direction(k, 0, -1, -1) + call s_add_particle_to_direction(k, -1, 0, -1) + call s_add_particle_to_direction(k, 0, 0, -1) + elseif (f_crosses_boundary(k, 3, 1, pos, posPrev)) then + call s_add_particle_to_direction(k, -1, -1, 1) + call s_add_particle_to_direction(k, 0, -1, 1) + call s_add_particle_to_direction(k, -1, 0, 1) + call s_add_particle_to_direction(k, 0, 0, 1) + end if + end if + elseif (f_crosses_boundary(k, 2, 1, pos, posPrev)) then + call s_add_particle_to_direction(k, -1, 1, 0) + call s_add_particle_to_direction(k, 0, 1, 0) + if (p > 0) then + if (f_crosses_boundary(k, 3, -1, pos, posPrev)) then + call s_add_particle_to_direction(k, -1, 1, -1) + call s_add_particle_to_direction(k, 0, 1, -1) + call s_add_particle_to_direction(k, -1, 0, -1) + call s_add_particle_to_direction(k, 0, 0, -1) + elseif (f_crosses_boundary(k, 3, 1, pos, posPrev)) then + call s_add_particle_to_direction(k, -1, 1, 1) + call s_add_particle_to_direction(k, 0, 1, 1) + call s_add_particle_to_direction(k, -1, 0, 1) + call s_add_particle_to_direction(k, 0, 0, 1) + end if + end if + else + if (p > 0) then + if (f_crosses_boundary(k, 3, -1, pos, posPrev)) then + call s_add_particle_to_direction(k, -1, 0, -1) + call s_add_particle_to_direction(k, 0, 0, -1) + elseif (f_crosses_boundary(k, 3, 1, pos, posPrev)) then + call s_add_particle_to_direction(k, -1, 0, 1) + call s_add_particle_to_direction(k, 0, 0, 1) + end if + end if + end if + end if + elseif (f_crosses_boundary(k, 1, 1, pos, posPrev)) then + call s_add_particle_to_direction(k, 1, 0, 0) + if (n > 0) then + if (f_crosses_boundary(k, 2, -1, pos, posPrev)) then + call s_add_particle_to_direction(k, 1, -1, 0) + call s_add_particle_to_direction(k, 0, -1, 0) + if (p > 0) then + if (f_crosses_boundary(k, 3, -1, pos, posPrev)) then + call s_add_particle_to_direction(k, 1, -1, -1) + call s_add_particle_to_direction(k, 0, -1, -1) + call s_add_particle_to_direction(k, 1, 0, -1) + call s_add_particle_to_direction(k, 0, 0, -1) + elseif (f_crosses_boundary(k, 3, 1, pos, posPrev)) then + call s_add_particle_to_direction(k, 1, -1, 1) + call s_add_particle_to_direction(k, 0, -1, 1) + call s_add_particle_to_direction(k, 1, 0, 1) + call s_add_particle_to_direction(k, 0, 0, 1) + end if + end if + elseif (f_crosses_boundary(k, 2, 1, pos, posPrev)) then + call s_add_particle_to_direction(k, 1, 1, 0) + call s_add_particle_to_direction(k, 0, 1, 0) + if (p > 0) then + if (f_crosses_boundary(k, 3, -1, pos, posPrev)) then + call s_add_particle_to_direction(k, 1, 1, -1) + call s_add_particle_to_direction(k, 0, 1, -1) + call s_add_particle_to_direction(k, 1, 0, -1) + call s_add_particle_to_direction(k, 0, 0, -1) + elseif (f_crosses_boundary(k, 3, 1, pos, posPrev)) then + call s_add_particle_to_direction(k, 1, 1, 1) + call s_add_particle_to_direction(k, 0, 1, 1) + call s_add_particle_to_direction(k, 1, 0, 1) + call s_add_particle_to_direction(k, 0, 0, 1) + end if + end if + else + if (p > 0) then + if (f_crosses_boundary(k, 3, -1, pos, posPrev)) then + call s_add_particle_to_direction(k, 1, 0, -1) + call s_add_particle_to_direction(k, 0, 0, -1) + elseif (f_crosses_boundary(k, 3, 1, pos, posPrev)) then + call s_add_particle_to_direction(k, 1, 0, 1) + call s_add_particle_to_direction(k, 0, 0, 1) + end if + end if + end if + end if + elseif (f_crosses_boundary(k, 2, -1, pos, posPrev)) then + call s_add_particle_to_direction(k, 0, -1, 0) + if (p > 0) then + if (f_crosses_boundary(k, 3, -1, pos, posPrev)) then + call s_add_particle_to_direction(k, 0, -1, -1) + call s_add_particle_to_direction(k, 0, 0, -1) + elseif (f_crosses_boundary(k, 3, 1, pos, posPrev)) then + call s_add_particle_to_direction(k, 0, -1, 1) + call s_add_particle_to_direction(k, 0, 0, 1) + end if + end if + elseif (f_crosses_boundary(k, 2, 1, pos, posPrev)) then + call s_add_particle_to_direction(k, 0, 1, 0) + if (p > 0) then + if (f_crosses_boundary(k, 3, -1, pos, posPrev)) then + call s_add_particle_to_direction(k, 0, 1, -1) + call s_add_particle_to_direction(k, 0, 0, -1) + elseif (f_crosses_boundary(k, 3, 1, pos, posPrev)) then + call s_add_particle_to_direction(k, 0, 1, 1) + call s_add_particle_to_direction(k, 0, 0, 1) + end if + end if + elseif (p > 0) then + if (f_crosses_boundary(k, 3, -1, pos, posPrev)) then + call s_add_particle_to_direction(k, 0, 0, -1) + elseif (f_crosses_boundary(k, 3, 1, pos, posPrev)) then + call s_add_particle_to_direction(k, 0, 0, 1) + end if + end if + + end do + + contains + + logical function f_crosses_boundary(particle_id, dir, loc, pos, posPrev) + + integer, intent(in) :: particle_id, dir, loc + real(wp), dimension(:, :), intent(in) :: pos + real(wp), dimension(:, :), optional, intent(in) :: posPrev + + if (loc == -1) then ! Beginning of the domain + if (nidx(dir)%beg == 0) then + f_crosses_boundary = .false. + return + end if + + if (present(posPrev)) then + f_crosses_boundary = (posPrev(particle_id, dir) > pcomm_coords(dir)%beg .and. & + pos(particle_id, dir) < pcomm_coords(dir)%beg) + else + f_crosses_boundary = (pos(particle_id, dir) < pcomm_coords(dir)%beg) + end if + elseif (loc == 1) then ! End of the domain + if (nidx(dir)%end == 0) then + f_crosses_boundary = .false. + return + end if + + if (present(posPrev)) then + f_crosses_boundary = (posPrev(particle_id, dir) < pcomm_coords(dir)%end .and. & + pos(particle_id, dir) > pcomm_coords(dir)%end) + else + f_crosses_boundary = (pos(particle_id, dir) > pcomm_coords(dir)%end) + end if + end if + + end function f_crosses_boundary + + subroutine s_add_particle_to_direction(particle_id, dir_x, dir_y, dir_z) + + integer, intent(in) :: particle_id, dir_x, dir_y, dir_z + + p_send_ids(dir_x, dir_y, dir_z, p_send_counts(dir_x, dir_y, dir_z)) = particle_id + p_send_counts(dir_x, dir_y, dir_z) = p_send_counts(dir_x, dir_y, dir_z) + 1 + + end subroutine s_add_particle_to_direction + + end subroutine s_add_particles_to_transfer_list + + !> This subroutine performs the MPI communication for lagrangian particles/ + !! bubbles. + !! @param bub_R0 Initial radius of each bubble + !! @param Rmax_stats Maximum radius of each bubble + !! @param Rmin_stats Minimum radius of each bubble + !! @param gas_mg Mass of gas in each bubble + !! @param gas_betaT Heat flux model coefficient for each bubble + !! @param gas_betaC mass flux model coefficient for each bubble + !! @param bub_dphidt Subgrid velocity potential for each bubble + !! @param lag_id Global and local ID of each bubble + !! @param gas_p Pressure of the gas in each bubble + !! @param gas_mv Mass of vapor in each bubble + !! @param rad Radius of each bubble + !! @param rvel Radial velocity of each bubble + !! @param pos Position of each bubble + !! @param posPrev Previous position of each bubble + !! @param vel Velocity of each bubble + !! @param scoord Cell index in real format of each bubble + !! @param drad Radial velocity of each bubble + !! @param drvel Radial acceleration of each bubble + !! @param dgasp Time derivative of gas pressure in each bubble + !! @param dgasmv Time derivative of vapor mass in each bubble + !! @param dpos Time derivative of position of each bubble + !! @param dvel Time derivative of velocity of each bubble + !! @param lag_num_ts Number of stages in time-stepping scheme + !! @param nBubs Local number of bubbles + impure subroutine s_mpi_sendrecv_particles(bub_R0, Rmax_stats, Rmin_stats, gas_mg, gas_betaT, & + gas_betaC, bub_dphidt, lag_id, gas_p, gas_mv, rad, & + rvel, pos, posPrev, vel, scoord, drad, drvel, dgasp, & + dgasmv, dpos, dvel, lag_num_ts, nbubs, dest) + + real(wp), dimension(:) :: bub_R0, Rmax_stats, Rmin_stats, gas_mg, gas_betaT, gas_betaC, bub_dphidt + integer, dimension(:, :) :: lag_id + real(wp), dimension(:, :) :: gas_p, gas_mv, rad, rvel, drad, drvel, dgasp, dgasmv + real(wp), dimension(:, :, :) :: pos, posPrev, vel, scoord, dpos, dvel + integer :: position, bub_id, lag_num_ts, tag, partner, send_tag, recv_tag, nbubs, p_recv_size, dest + + integer :: i, j, k, l, q, r + integer :: req_send, req_recv, ierr !< Generic flag used to identify and report MPI errors + integer :: send_count, send_offset, recv_count, recv_offset + +#ifdef MFC_MPI + ! Phase 1: Exchange particle counts using non-blocking communication + send_count = 0 + recv_count = 0 + + ! Post all receives first + do l = 1, n_neighbors + i = neighbor_list(l, 1) + j = neighbor_list(l, 2) + k = neighbor_list(l, 3) + partner = neighbor_ranks(i, j, k) + recv_tag = neighbor_tag(i, j, k) + + recv_count = recv_count + 1 + call MPI_Irecv(p_recv_counts(i, j, k), 1, MPI_INTEGER, partner, recv_tag, & + MPI_COMM_WORLD, recv_requests(recv_count), ierr) + end do + + ! Post all sends + do l = 1, n_neighbors + i = neighbor_list(l, 1) + j = neighbor_list(l, 2) + k = neighbor_list(l, 3) + partner = neighbor_ranks(i, j, k) + send_tag = neighbor_tag(-i, -j, -k) + + send_count = send_count + 1 + call MPI_Isend(p_send_counts(i, j, k), 1, MPI_INTEGER, partner, send_tag, & + MPI_COMM_WORLD, send_requests(send_count), ierr) + end do + + ! Wait for all count exchanges to complete + if (recv_count > 0) then + call MPI_Waitall(recv_count, recv_requests(1:recv_count), MPI_STATUSES_IGNORE, ierr) + end if + if (send_count > 0) then + call MPI_Waitall(send_count, send_requests(1:send_count), MPI_STATUSES_IGNORE, ierr) + end if + + ! Phase 2: Exchange particle data using non-blocking communication + send_count = 0 + recv_count = 0 + + ! Post all receives for particle data first + recv_offset = 1 + do l = 1, n_neighbors + i = neighbor_list(l, 1) + j = neighbor_list(l, 2) + k = neighbor_list(l, 3) + + if (p_recv_counts(i, j, k) > 0) then + partner = neighbor_ranks(i, j, k) + p_recv_size = p_recv_counts(i, j, k) * p_var_size + recv_tag = neighbor_tag(i, j, k) + + recv_count = recv_count + 1 + call MPI_Irecv(p_recv_buff(recv_offset), p_recv_size, MPI_PACKED, partner, recv_tag, & + MPI_COMM_WORLD, recv_requests(recv_count), ierr) + recv_offsets(l) = recv_offset + recv_offset = recv_offset + p_recv_size + end if + end do + + ! Pack and send particle data + send_offset = 0 + do l = 1, n_neighbors + i = neighbor_list(l, 1) + j = neighbor_list(l, 2) + k = neighbor_list(l, 3) + + if (p_send_counts(i, j, k) > 0) then + partner = neighbor_ranks(i, j, k) + send_tag = neighbor_tag(-i, -j, -k) + + ! Pack data for sending + position = 0 + do q = 0, p_send_counts(i, j, k) - 1 + bub_id = p_send_ids(i, j, k, q) + call MPI_Pack(lag_id(bub_id, 1), 1, MPI_INTEGER, p_send_buff(send_offset), p_buff_size, position, MPI_COMM_WORLD, ierr) + call MPI_Pack(bub_R0(bub_id), 1, mpi_p, p_send_buff(send_offset), p_buff_size, position, MPI_COMM_WORLD, ierr) + call MPI_Pack(Rmax_stats(bub_id), 1, mpi_p, p_send_buff(send_offset), p_buff_size, position, MPI_COMM_WORLD, ierr) + call MPI_Pack(Rmin_stats(bub_id), 1, mpi_p, p_send_buff(send_offset), p_buff_size, position, MPI_COMM_WORLD, ierr) + call MPI_Pack(gas_mg(bub_id), 1, mpi_p, p_send_buff(send_offset), p_buff_size, position, MPI_COMM_WORLD, ierr) + call MPI_Pack(gas_betaT(bub_id), 1, mpi_p, p_send_buff(send_offset), p_buff_size, position, MPI_COMM_WORLD, ierr) + call MPI_Pack(gas_betaC(bub_id), 1, mpi_p, p_send_buff(send_offset), p_buff_size, position, MPI_COMM_WORLD, ierr) + call MPI_Pack(bub_dphidt(bub_id), 1, mpi_p, p_send_buff(send_offset), p_buff_size, position, MPI_COMM_WORLD, ierr) + do r = 1, 2 + call MPI_Pack(gas_p(bub_id, r), 1, mpi_p, p_send_buff(send_offset), p_buff_size, position, MPI_COMM_WORLD, ierr) + call MPI_Pack(gas_mv(bub_id, r), 1, mpi_p, p_send_buff(send_offset), p_buff_size, position, MPI_COMM_WORLD, ierr) + call MPI_Pack(rad(bub_id, r), 1, mpi_p, p_send_buff(send_offset), p_buff_size, position, MPI_COMM_WORLD, ierr) + call MPI_Pack(rvel(bub_id, r), 1, mpi_p, p_send_buff(send_offset), p_buff_size, position, MPI_COMM_WORLD, ierr) + call MPI_Pack(pos(bub_id, :, r), 3, mpi_p, p_send_buff(send_offset), p_buff_size, position, MPI_COMM_WORLD, ierr) + call MPI_Pack(posPrev(bub_id, :, r), 3, mpi_p, p_send_buff(send_offset), p_buff_size, position, MPI_COMM_WORLD, ierr) + call MPI_Pack(vel(bub_id, :, r), 3, mpi_p, p_send_buff(send_offset), p_buff_size, position, MPI_COMM_WORLD, ierr) + call MPI_Pack(scoord(bub_id, :, r), 3, mpi_p, p_send_buff(send_offset), p_buff_size, position, MPI_COMM_WORLD, ierr) + end do + do r = 1, lag_num_ts + call MPI_Pack(drad(bub_id, r), 1, mpi_p, p_send_buff(send_offset), p_buff_size, position, MPI_COMM_WORLD, ierr) + call MPI_Pack(drvel(bub_id, r), 1, mpi_p, p_send_buff(send_offset), p_buff_size, position, MPI_COMM_WORLD, ierr) + call MPI_Pack(dgasp(bub_id, r), 1, mpi_p, p_send_buff(send_offset), p_buff_size, position, MPI_COMM_WORLD, ierr) + call MPI_Pack(dgasmv(bub_id, r), 1, mpi_p, p_send_buff(send_offset), p_buff_size, position, MPI_COMM_WORLD, ierr) + call MPI_Pack(dpos(bub_id, :, r), 3, mpi_p, p_send_buff(send_offset), p_buff_size, position, MPI_COMM_WORLD, ierr) + call MPI_Pack(dvel(bub_id, :, r), 3, mpi_p, p_send_buff(send_offset), p_buff_size, position, MPI_COMM_WORLD, ierr) + end do + end do + + send_count = send_count + 1 + call MPI_Isend(p_send_buff(send_offset), position, MPI_PACKED, partner, send_tag, & + MPI_COMM_WORLD, send_requests(send_count), ierr) + send_offset = send_offset + position + end if + end do + + ! Wait for all recvs for contiguous data to complete + call MPI_Waitall(recv_count, recv_requests(1:recv_count), MPI_STATUSES_IGNORE, ierr) + + ! Process received data as it arrives + do l = 1, n_neighbors + i = neighbor_list(l, 1) + j = neighbor_list(l, 2) + k = neighbor_list(l, 3) + + if (p_recv_counts(i, j, k) > 0 .and. abs(i) + abs(j) + abs(k) /= 0) then + p_recv_size = p_recv_counts(i, j, k) * p_var_size + recv_offset = recv_offsets(l) + + position = 0 + ! Unpack received data + do q = 0, p_recv_counts(i, j, k) - 1 + nbubs = nbubs + 1 + bub_id = nbubs + call MPI_Unpack(p_recv_buff(recv_offset), p_recv_size, position, lag_id(bub_id, 1), 1, MPI_INTEGER, MPI_COMM_WORLD, ierr) + call MPI_Unpack(p_recv_buff(recv_offset), p_recv_size, position, bub_R0(bub_id), 1, mpi_p, MPI_COMM_WORLD, ierr) + call MPI_Unpack(p_recv_buff(recv_offset), p_recv_size, position, Rmax_stats(bub_id), 1, mpi_p, MPI_COMM_WORLD, ierr) + call MPI_Unpack(p_recv_buff(recv_offset), p_recv_size, position, Rmin_stats(bub_id), 1, mpi_p, MPI_COMM_WORLD, ierr) + call MPI_Unpack(p_recv_buff(recv_offset), p_recv_size, position, gas_mg(bub_id), 1, mpi_p, MPI_COMM_WORLD, ierr) + call MPI_Unpack(p_recv_buff(recv_offset), p_recv_size, position, gas_betaT(bub_id), 1, mpi_p, MPI_COMM_WORLD, ierr) + call MPI_Unpack(p_recv_buff(recv_offset), p_recv_size, position, gas_betaC(bub_id), 1, mpi_p, MPI_COMM_WORLD, ierr) + call MPI_Unpack(p_recv_buff(recv_offset), p_recv_size, position, bub_dphidt(bub_id), 1, mpi_p, MPI_COMM_WORLD, ierr) + do r = 1, 2 + call MPI_Unpack(p_recv_buff(recv_offset), p_recv_size, position, gas_p(bub_id, r), 1, mpi_p, MPI_COMM_WORLD, ierr) + call MPI_Unpack(p_recv_buff(recv_offset), p_recv_size, position, gas_mv(bub_id, r), 1, mpi_p, MPI_COMM_WORLD, ierr) + call MPI_Unpack(p_recv_buff(recv_offset), p_recv_size, position, rad(bub_id, r), 1, mpi_p, MPI_COMM_WORLD, ierr) + call MPI_Unpack(p_recv_buff(recv_offset), p_recv_size, position, rvel(bub_id, r), 1, mpi_p, MPI_COMM_WORLD, ierr) + call MPI_Unpack(p_recv_buff(recv_offset), p_recv_size, position, pos(bub_id, :, r), 3, mpi_p, MPI_COMM_WORLD, ierr) + call MPI_Unpack(p_recv_buff(recv_offset), p_recv_size, position, posPrev(bub_id, :, r), 3, mpi_p, MPI_COMM_WORLD, ierr) + call MPI_Unpack(p_recv_buff(recv_offset), p_recv_size, position, vel(bub_id, :, r), 3, mpi_p, MPI_COMM_WORLD, ierr) + call MPI_Unpack(p_recv_buff(recv_offset), p_recv_size, position, scoord(bub_id, :, r), 3, mpi_p, MPI_COMM_WORLD, ierr) + end do + do r = 1, lag_num_ts + call MPI_Unpack(p_recv_buff(recv_offset), p_recv_size, position, drad(bub_id, r), 1, mpi_p, MPI_COMM_WORLD, ierr) + call MPI_Unpack(p_recv_buff(recv_offset), p_recv_size, position, drvel(bub_id, r), 1, mpi_p, MPI_COMM_WORLD, ierr) + call MPI_Unpack(p_recv_buff(recv_offset), p_recv_size, position, dgasp(bub_id, r), 1, mpi_p, MPI_COMM_WORLD, ierr) + call MPI_Unpack(p_recv_buff(recv_offset), p_recv_size, position, dgasmv(bub_id, r), 1, mpi_p, MPI_COMM_WORLD, ierr) + call MPI_Unpack(p_recv_buff(recv_offset), p_recv_size, position, dpos(bub_id, :, r), 3, mpi_p, MPI_COMM_WORLD, ierr) + call MPI_Unpack(p_recv_buff(recv_offset), p_recv_size, position, dvel(bub_id, :, r), 3, mpi_p, MPI_COMM_WORLD, ierr) + end do + lag_id(bub_id, 2) = bub_id + end do + recv_offset = recv_offset + p_recv_size + end if + + end do + + ! Wait for all sends to complete + if (send_count > 0) then + call MPI_Waitall(send_count, send_requests(1:send_count), MPI_STATUSES_IGNORE, ierr) + end if +#endif + + if (any(periodic_bc)) then + call s_wrap_particle_positions(pos, posPrev, nbubs, dest) + end if + + end subroutine s_mpi_sendrecv_particles + + !! This function returns a unique tag for each neighbor based on its position + !! relative to the current process. + !! @param i, j, k Indices of the neighbor in the range [-1, 1] + !! @return tag Unique integer tag for the neighbor + integer function neighbor_tag(i, j, k) result(tag) + + integer, intent(in) :: i, j, k + + tag = (k + 1)*9 + (j + 1)*3 + (i + 1) + + end function neighbor_tag + + subroutine s_wrap_particle_positions(pos, posPrev, nbubs, dest) + + real(wp), dimension(:, :, :) :: pos, posPrev + integer :: nbubs, dest + integer :: i, q + real :: offset + + do i = 1, nbubs + if (periodic_bc(1)) then + offset = glb_bounds(1)%end - glb_bounds(1)%beg + if (pos(i, 1, dest) > x_cb(m + buff_size)) then + do q = 1, 2 + pos(i, 1, q) = pos(i, 1, q) - offset + posPrev(i, 1, q) = posPrev(i, 1, q) - offset + end do + endif + if (pos(i, 1, dest) < x_cb(-1 - buff_size)) then + do q = 1, 2 + pos(i, 1, q) = pos(i, 1, q) + offset + posPrev(i, 1, q) = posPrev(i, 1, q) + offset + end do + endif + end if + + if (periodic_bc(2)) then + offset = glb_bounds(2)%end - glb_bounds(2)%beg + if (pos(i, 2, dest) > y_cb(n + buff_size)) then + do q = 1, 2 + pos(i, 2, q) = pos(i, 2, q) - offset + posPrev(i, 2, q) = posPrev(i, 2, q) - offset + end do + endif + if (pos(i, 2, dest) < y_cb(-buff_size - 1)) then + do q = 1, 2 + pos(i, 2, q) = pos(i, 2, q) + offset + posPrev(i, 2, q) = posPrev(i, 2, q) + offset + end do + endif + end if + + if (periodic_bc(3)) then + offset = glb_bounds(3)%end - glb_bounds(3)%beg + if (pos(i, 3, dest) > z_cb(p + buff_size)) then + do q = 1, 2 + pos(i, 3, q) = pos(i, 3, q) - offset + posPrev(i, 3, q) = posPrev(i, 3, q) - offset + end do + endif + if (pos(i, 3, dest) < z_cb(-1 - buff_size)) then + do q = 1, 2 + pos(i, 3, q) = pos(i, 3, q) + offset + posPrev(i, 3, q) = posPrev(i, 3, q) + offset + end do + endif + end if + end do + + end subroutine s_wrap_particle_positions + impure subroutine s_mpi_send_random_number(phi_rn, num_freq) integer, intent(in) :: num_freq real(wp), intent(inout), dimension(1:num_freq) :: phi_rn diff --git a/src/simulation/m_rhs.fpp b/src/simulation/m_rhs.fpp index cdcb418075..3a17174b72 100644 --- a/src/simulation/m_rhs.fpp +++ b/src/simulation/m_rhs.fpp @@ -1003,13 +1003,6 @@ contains end if if (bubbles_lagrange) then - ! RHS additions for sub-grid bubbles_lagrange - call nvtxStartRange("RHS-EL-BUBBLES-SRC") - call s_compute_bubbles_EL_source( & - q_cons_qp%vf(1:sys_size), & - q_prim_qp%vf(1:sys_size), & - rhs_vf) - call nvtxEndRange ! Compute bubble dynamics if (.not. adap_dt) then call nvtxStartRange("RHS-EL-BUBBLES-DYN") @@ -1018,6 +1011,14 @@ contains stage) call nvtxEndRange end if + + ! RHS additions for sub-grid bubbles_lagrange + call nvtxStartRange("RHS-EL-BUBBLES-SRC") + call s_compute_bubbles_EL_source( & + q_cons_qp%vf(1:sys_size), & + q_prim_qp%vf(1:sys_size), & + rhs_vf) + call nvtxEndRange end if if (chemistry .and. chem_params%reactions) then diff --git a/src/simulation/m_sim_helpers.fpp b/src/simulation/m_sim_helpers.fpp index df8f6b9eae..01fc7a9ec1 100644 --- a/src/simulation/m_sim_helpers.fpp +++ b/src/simulation/m_sim_helpers.fpp @@ -181,12 +181,12 @@ contains !! @param icfl_sf cell-centered inviscid cfl number !! @param vcfl_sf (optional) cell-centered viscous CFL number !! @param Rc_sf (optional) cell centered Rc - pure subroutine s_compute_stability_from_dt(vel, c, rho, Re_l, j, k, l, icfl_sf, vcfl_sf, Rc_sf) + pure subroutine s_compute_stability_from_dt(vel, c, rho, Re_l, j, k, l, icfl_sf, vcfl_sf, Rc_sf, ccfl_sf) $:GPU_ROUTINE(parallelism='[seq]') real(wp), intent(in), dimension(num_vels) :: vel real(wp), intent(in) :: c, rho real(wp), dimension(0:m, 0:n, 0:p), intent(inout) :: icfl_sf - real(wp), dimension(0:m, 0:n, 0:p), intent(inout), optional :: vcfl_sf, Rc_sf + real(wp), dimension(0:m, 0:n, 0:p), intent(inout), optional :: vcfl_sf, Rc_sf, ccfl_sf real(wp), dimension(2), intent(in) :: Re_l integer, intent(in) :: j, k, l @@ -203,8 +203,7 @@ contains ! Viscous calculations if (viscous) then - if (p > 0) then - !3D + if (p > 0) then ! 3D if (grid_geometry == 3) then fltr_dtheta = f_compute_filtered_dtheta(k, l) vcfl_sf(j, k, l) = maxval(dt/Re_l/rho) & @@ -221,19 +220,30 @@ contains dz(l)*(abs(vel(3)) + c)) & /maxval(1._wp/Re_l) end if - elseif (n > 0) then - !2D + elseif (n > 0) then ! 2D vcfl_sf(j, k, l) = maxval(dt/Re_l/rho)/min(dx(j), dy(k))**2._wp Rc_sf(j, k, l) = min(dx(j)*(abs(vel(1)) + c), & dy(k)*(abs(vel(2)) + c)) & /maxval(1._wp/Re_l) - else - !1D + else ! 1D vcfl_sf(j, k, l) = maxval(dt/Re_l/rho)/dx(j)**2._wp Rc_sf(j, k, l) = dx(j)*(abs(vel(1)) + c)/maxval(1._wp/Re_l) end if end if + if (surface_tension) then + if (p > 0) then ! 3D + if (grid_geometry == 3) then + fltr_dtheta = f_compute_filtered_dtheta(k, l) + ccfl_sf(j, k, l) = dt/sqrt(rho * min(dx(j), dy(k), fltr_dtheta)**3._wp / sigma) + else + ccfl_sf(j, k, l) = dt/sqrt(rho * min(dx(j), dy(k), dz(l))**3._wp / sigma) + end if + elseif (n > 0) then ! 2D + ccfl_sf(j, k, l) = dt/sqrt(rho * min(dx(j), dy(k))**3._wp / sigma) + end if + end if + end subroutine s_compute_stability_from_dt !> Computes dt for a specified CFL number @@ -252,7 +262,7 @@ contains real(wp), dimension(2), intent(in) :: Re_l integer, intent(in) :: j, k, l - real(wp) :: icfl_dt, vcfl_dt + real(wp) :: icfl_dt, vcfl_dt, ccfl_dt real(wp) :: fltr_dtheta ! Inviscid CFL calculation @@ -266,8 +276,7 @@ contains ! Viscous calculations if (viscous) then - if (p > 0) then - !3D + if (p > 0) then ! 3D if (grid_geometry == 3) then fltr_dtheta = f_compute_filtered_dtheta(k, l) vcfl_dt = cfl_target*(min(dx(j), dy(k), fltr_dtheta)**2._wp) & @@ -276,17 +285,32 @@ contains vcfl_dt = cfl_target*(min(dx(j), dy(k), dz(l))**2._wp) & /minval(1/(rho*Re_l)) end if - elseif (n > 0) then - !2D + elseif (n > 0) then ! 2D vcfl_dt = cfl_target*(min(dx(j), dy(k))**2._wp)/maxval((1/Re_l)/rho) - else - !1D + else ! 1D vcfl_dt = cfl_target*(dx(j)**2._wp)/minval(1/(rho*Re_l)) end if end if - if (any(Re_size > 0)) then + if (surface_tension) then + if (p > 0) then ! 3D + if (grid_geometry == 3) then + fltr_dtheta = f_compute_filtered_dtheta(k, l) + ccfl_dt = cfl_target*sqrt(rho * min(dx(j), dy(k), fltr_dtheta)**3._wp / sigma) + else + ccfl_dt = cfl_target*sqrt(rho * min(dx(j), dy(k), dz(l))**3._wp / sigma) + end if + elseif (n > 0) then ! 2D + ccfl_dt = cfl_target*sqrt(rho * min(dx(j), dy(k))**3._wp / sigma) + end if + end if + + if (any(Re_size > 0) .and. sigma > 0) then + max_dt(j, k, l) = min(icfl_dt, vcfl_dt, ccfl_dt) + elseif (any(Re_size > 0)) then max_dt(j, k, l) = min(icfl_dt, vcfl_dt) + elseif (sigma > 0) then + max_dt(j, k, l) = min(icfl_dt, ccfl_dt) else max_dt(j, k, l) = icfl_dt end if diff --git a/src/simulation/m_start_up.fpp b/src/simulation/m_start_up.fpp index 6705c68bf4..da4bb870b4 100644 --- a/src/simulation/m_start_up.fpp +++ b/src/simulation/m_start_up.fpp @@ -155,8 +155,6 @@ contains rdma_mpi, teno_CT, mp_weno, weno_avg, & riemann_solver, low_Mach, wave_speeds, avg_state, & bc_x, bc_y, bc_z, & - x_a, y_a, z_a, x_b, y_b, z_b, & - x_domain, y_domain, z_domain, & hypoelasticity, & ib, num_ibs, patch_ib, & fluid_pp, probe_wrt, prim_vars_wrt, & @@ -225,11 +223,15 @@ contains if (cfl_adap_dt .or. cfl_const_dt) cfl_dt = .true. - if (any((/bc_x%beg, bc_x%end, bc_y%beg, bc_y%end, bc_z%beg, bc_z%end/) == -17) .or. & + if (any((/bc_x%beg, bc_x%end, bc_y%beg, bc_y%end, bc_z%beg, bc_z%end/) == BC_DIRICHLET) .or. & num_bc_patches > 0) then bc_io = .true. endif + if (bc_x%beg == BC_PERIODIC .and. bc_x%end == BC_PERIODIC) periodic_bc(1) = .true. + if (bc_y%beg == BC_PERIODIC .and. bc_y%end == BC_PERIODIC) periodic_bc(2) = .true. + if (bc_z%beg == BC_PERIODIC .and. bc_z%end == BC_PERIODIC) periodic_bc(3) = .true. + else call s_mpi_abort(trim(file_path)//' is missing. Exiting.') end if @@ -1252,7 +1254,9 @@ contains end if if (bubbles_lagrange) then - $:GPU_UPDATE(host='[intfc_rad]') + $:GPU_UPDATE(host='[lag_id, mtn_pos, mtn_posPrev, mtn_vel, intfc_rad, & + & intfc_vel, bub_R0, Rmax_stats, Rmin_stats, bub_dphidt, gas_p, & + & gas_mv, gas_mg, gas_betaT, gas_betaC]') do i = 1, nBubs if (ieee_is_nan(intfc_rad(i, 1)) .or. intfc_rad(i, 1) <= 0._wp) then call s_mpi_abort("Bubble radius is negative or NaN, please reduce dt.") @@ -1261,8 +1265,9 @@ contains $:GPU_UPDATE(host='[q_beta%vf(1)%sf]') call s_write_data_files(q_cons_ts(1)%vf, q_T_sf, q_prim_vf, save_count, bc_type, q_beta%vf(1)) - $:GPU_UPDATE(host='[Rmax_stats,Rmin_stats,gas_p,gas_mv,intfc_vel]') + call s_write_restart_lag_bubbles(save_count) !parallel + if (lag_params%write_bubbles_stats) call s_write_lag_bubble_stats() else call s_write_data_files(q_cons_ts(1)%vf, q_T_sf, q_prim_vf, save_count, bc_type) @@ -1384,8 +1389,8 @@ contains end if call s_initialize_derived_variables() - if (bubbles_lagrange) call s_initialize_bubbles_EL_module(q_cons_ts(1)%vf) + if (bubbles_lagrange) call s_initialize_bubbles_EL_module(q_cons_ts(1)%vf) if (hypoelasticity) call s_initialize_hypoelastic_module() if (hyperelasticity) call s_initialize_hyperelastic_module() @@ -1501,6 +1506,10 @@ contains $:GPU_UPDATE(device='[bc_y%vb1,bc_y%vb2,bc_y%vb3,bc_y%ve1,bc_y%ve2,bc_y%ve3]') $:GPU_UPDATE(device='[bc_z%vb1,bc_z%vb2,bc_z%vb3,bc_z%ve1,bc_z%ve2,bc_z%ve3]') + $:GPU_UPDATE(device='[bc_x%beg, bc_x%end]') + $:GPU_UPDATE(device='[bc_y%beg, bc_y%end]') + $:GPU_UPDATE(device='[bc_z%beg, bc_z%end]') + $:GPU_UPDATE(device='[bc_x%grcbc_in,bc_x%grcbc_out,bc_x%grcbc_vel_out]') $:GPU_UPDATE(device='[bc_y%grcbc_in,bc_y%grcbc_out,bc_y%grcbc_vel_out]') $:GPU_UPDATE(device='[bc_z%grcbc_in,bc_z%grcbc_out,bc_z%grcbc_vel_out]') diff --git a/src/simulation/m_time_steppers.fpp b/src/simulation/m_time_steppers.fpp index d272c842aa..ae4fe4aeb4 100644 --- a/src/simulation/m_time_steppers.fpp +++ b/src/simulation/m_time_steppers.fpp @@ -506,7 +506,7 @@ contains end if end if - if (bubbles_lagrange .and. .not. adap_dt) call s_update_lagrange_tdv_rk(stage=s) + if (bubbles_lagrange .and. .not. adap_dt) call s_update_lagrange_tdv_rk(q_prim_vf, stage=s) $:GPU_PARALLEL_LOOP(collapse=4) do i = 1, sys_size @@ -616,7 +616,6 @@ contains else wall_time_avg = 0._wp end if - end subroutine s_tvd_rk !> Bubble source part in Strang operator splitting scheme diff --git a/toolchain/mfc/run/case_dicts.py b/toolchain/mfc/run/case_dicts.py index 044e933f3f..9609407c78 100644 --- a/toolchain/mfc/run/case_dicts.py +++ b/toolchain/mfc/run/case_dicts.py @@ -106,6 +106,13 @@ def analytic(self): 'elliptic_smoothing_iters': ParamType.INT, 'viscous': ParamType.LOG, 'bubbles_lagrange': ParamType.LOG, + 'lag_params%vel_model': ParamType.INT, + 'fd_order': ParamType.INT, + 'normFac': ParamType.REAL, + 'interface_file': ParamType.STR, + 'normMag': ParamType.REAL, + 'g0': ParamType.REAL, + 'p0': ParamType.REAL, }) for ib_id in range(1, 10+1): @@ -162,7 +169,7 @@ def analytic(self): PRE_PROCESS[f"patch_bc({bc_p_id})%radius"] = ParamType.REAL -for p_id in range(1, 10+1): +for p_id in range(1, 20+1): for attribute, ty in [("geometry", ParamType.INT), ("smoothen", ParamType.LOG), ("smooth_patch_id", ParamType.INT), ("hcid", ParamType.INT)]: PRE_PROCESS[f"patch_icpp({p_id})%{attribute}"] = ty @@ -325,14 +332,16 @@ def analytic(self): }) for var in [ 'heatTransfer_model', 'massTransfer_model', 'pressure_corrector', - 'write_bubbles', 'write_bubbles_stats' ]: + 'write_bubbles', 'write_bubbles_stats', 'pressure_force', + 'gravity_force', 'momentum_transfer_force']: SIMULATION[f'lag_params%{var}'] = ParamType.LOG -for var in [ 'solver_approach', 'cluster_type', 'smooth_type', 'nBubs_glb']: +for var in [ 'solver_approach', 'cluster_type', 'smooth_type', 'nBubs_glb', + 'vel_model', 'drag_model']: SIMULATION[f'lag_params%{var}'] = ParamType.INT for var in [ 'epsilonb', 'valmaxvoid', 'charwidth', 'diffcoefvap', - 'c0', 'rho0', 'T0', 'x0', 'Thost' ]: + 'c0', 'rho0', 'T0', 'x0', 'Thost', 'c_d' ]: SIMULATION[f'lag_params%{var}'] = ParamType.REAL for var in [ 'diffusion', 'reactions' ]: @@ -386,10 +395,6 @@ def analytic(self): SIMULATION[f'{var}_{cmp}'] = ParamType.REAL SIMULATION[f'bf_{cmp}'] = ParamType.LOG - - for prepend in ["domain%beg", "domain%end"]: - SIMULATION[f"{cmp}_{prepend}"] = ParamType.REAL - for probe_id in range(1,10+1): for cmp in ["x", "y", "z"]: SIMULATION[f'probe({probe_id})%{cmp}'] = ParamType.REAL @@ -452,6 +457,24 @@ def analytic(self): 'pres_inf_wrt': ParamType.LOG, 'cons_vars_wrt': ParamType.LOG, 'prim_vars_wrt': ParamType.LOG, + 'lag_header': ParamType.LOG, + 'lag_txt_wrt': ParamType.LOG, + 'lag_db_wrt': ParamType.LOG, + 'lag_id_wrt': ParamType.LOG, + 'lag_pos_wrt': ParamType.LOG, + 'lag_pos_prev_wrt': ParamType.LOG, + 'lag_vel_wrt': ParamType.LOG, + 'lag_rad_wrt': ParamType.LOG, + 'lag_rvel_wrt': ParamType.LOG, + 'lag_r0_wrt': ParamType.LOG, + 'lag_rmax_wrt': ParamType.LOG, + 'lag_rmin_wrt': ParamType.LOG, + 'lag_dphidt_wrt': ParamType.LOG, + 'lag_pres_wrt': ParamType.LOG, + 'lag_mv_wrt': ParamType.LOG, + 'lag_mg_wrt': ParamType.LOG, + 'lag_betaT_wrt': ParamType.LOG, + 'lag_betaC_wrt': ParamType.LOG, 'c_wrt': ParamType.LOG, 'omega_wrt': ParamType.LOG, 'qbmm': ParamType.LOG, @@ -468,6 +491,24 @@ def analytic(self): 'surface_tension': ParamType.LOG, 'output_partial_domain': ParamType.LOG, 'bubbles_lagrange': ParamType.LOG, + 'lag_header': ParamType.LOG, + 'lag_txt_wrt': ParamType.LOG, + 'lag_db_wrt': ParamType.LOG, + 'lag_id_wrt': ParamType.LOG, + 'lag_pos_wrt': ParamType.LOG, + 'lag_pos_prev_wrt': ParamType.LOG, + 'lag_vel_wrt': ParamType.LOG, + 'lag_rad_wrt': ParamType.LOG, + 'lag_rvel_wrt': ParamType.LOG, + 'lag_r0_wrt': ParamType.LOG, + 'lag_rmax_wrt': ParamType.LOG, + 'lag_rmin_wrt': ParamType.LOG, + 'lag_dphidt_wrt': ParamType.LOG, + 'lag_pres_wrt': ParamType.LOG, + 'lag_mv_wrt': ParamType.LOG, + 'lag_mg_wrt': ParamType.LOG, + 'lag_betaT_wrt': ParamType.LOG, + 'lag_betaC_wrt': ParamType.LOG, }) for cmp in ["x", "y", "z"]: