diff --git a/.github/workflows/test.yml b/.github/workflows/test.yml index 3031858eb1..503f691ec6 100644 --- a/.github/workflows/test.yml +++ b/.github/workflows/test.yml @@ -48,6 +48,11 @@ jobs: - name: Clone uses: actions/checkout@v4 + - name: Set up Python 3.13 + uses: actions/setup-python@v5 + with: + python-version: '3.13' + - name: Setup MacOS if: matrix.os == 'macos' run: | 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/m_constants.fpp b/src/common/m_constants.fpp index 1d2e53d206..114286f53b 100644 --- a/src/common/m_constants.fpp +++ b/src/common/m_constants.fpp @@ -61,6 +61,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_mpi_common.fpp b/src/common/m_mpi_common.fpp index 07463137f2..bc63e1bdfc 100644 --- a/src/common/m_mpi_common.fpp +++ b/src/common/m_mpi_common.fpp @@ -1570,14 +1570,14 @@ contains #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 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..07046fc2bf 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 @@ -258,6 +257,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 +346,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 @@ -438,6 +457,24 @@ contains sim_data = .false. cf_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 +849,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 +858,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 +899,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..7ef6de571d 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 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/simulation/m_body_forces.fpp b/src/simulation/m_body_forces.fpp index 41a103fcee..3b88efe060 100644 --- a/src/simulation/m_body_forces.fpp +++ b/src/simulation/m_body_forces.fpp @@ -57,15 +57,11 @@ contains real(wp), intent(in) :: t - 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..23c255428a 100644 --- a/src/simulation/m_bubbles.fpp +++ b/src/simulation/m_bubbles.fpp @@ -72,6 +72,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 +448,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 diff --git a/src/simulation/m_bubbles_EL.fpp b/src/simulation/m_bubbles_EL.fpp index f5dbc1c628..7d0dd32bd8 100644 --- a/src/simulation/m_bubbles_EL.fpp +++ b/src/simulation/m_bubbles_EL.fpp @@ -226,7 +226,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") @@ -324,25 +324,21 @@ contains cell = -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 +401,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 +417,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)) + + call MPI_TYPE_CREATE_SUBARRAY(2, gsizes, lsizes, start_idx_part, & + MPI_ORDER_FORTRAN, mpi_p, view, ierr) + call MPI_TYPE_COMMIT(view, ierr) - if (particle_file) then - call MPI_FILE_open(MPI_COMM_WORLD, file_loc, MPI_MODE_RDONLY, & + 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 + + 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 - call MPI_FILE_CLOSE(ifile, ierr) + + deallocate (proc_bubble_counts) #endif end subroutine s_restart_bubbles @@ -649,6 +712,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 +750,7 @@ contains call s_gradient_dir(q_prim_vf(E_idx), q_beta%vf(3), l) + ! (q / (1 - beta)) * d(beta)/dt source $:GPU_PARALLEL_LOOP(collapse=3) do k = 0, p do j = 0, n @@ -712,6 +777,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 @@ -1234,25 +1300,24 @@ contains (pos_part(3) < z_cb(p + buff_size)) .and. (pos_part(3) >= z_cb(-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 +1368,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 +1528,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 +1539,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 +1554,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_ALLGATHER(bub_id, 1, MPI_INTEGER, proc_bubble_counts, 1, MPI_INTEGER, & + MPI_COMM_WORLD, ierr) - call MPI_ALLREDUCE(part_order, part_ord_mpi, num_procs, MPI_integer, & - MPI_MAX, 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) - - allocate (MPI_IO_DATA_lag_bubbles(1:max(1, bub_id), 1:21)) + call MPI_BARRIER(MPI_COMM_WORLD, ierr) - ! 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 - - call MPI_FILE_OPEN(MPI_COMM_WORLD, file_loc, ior(MPI_MODE_WRONLY, MPI_MODE_CREATE), & - mpi_info_int, ifile, ierr) + 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) - disp = 0._wp + ! 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) - call MPI_FILE_SET_VIEW(ifile, disp, mpi_p, view, & - 'native', mpi_info_null, ierr) + call MPI_FILE_CLOSE(ifile, ierr) + end if - ! Cycle through list - i = 1 + call MPI_BARRIER(MPI_COMM_WORLD, ierr) - 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 +1623,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_CLOSE(ifile, ierr) - call MPI_FILE_write_ALL(ifile, MPI_IO_DATA_lag_bubbles, 21*max(1, bub_id), & - mpi_p, status, 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_CLOSE(ifile, 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 diff --git a/src/simulation/m_data_output.fpp b/src/simulation/m_data_output.fpp index 5f5b79c481..0bc804b183 100644 --- a/src/simulation/m_data_output.fpp +++ b/src/simulation/m_data_output.fpp @@ -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 @@ -358,16 +350,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 diff --git a/src/simulation/m_start_up.fpp b/src/simulation/m_start_up.fpp index 123db558ea..c83d22e697 100644 --- a/src/simulation/m_start_up.fpp +++ b/src/simulation/m_start_up.fpp @@ -1265,7 +1265,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.") @@ -1274,7 +1276,7 @@ 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 diff --git a/toolchain/mfc/run/case_dicts.py b/toolchain/mfc/run/case_dicts.py index 56eb21992c..b9c28717b1 100644 --- a/toolchain/mfc/run/case_dicts.py +++ b/toolchain/mfc/run/case_dicts.py @@ -464,6 +464,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"]: