Merge branch 'master' into raw-directives

Author: sbryngelson

README.md

@@ -137,7 +137,6 @@ They are organized below.
 * Acoustic wave generation (one- and two-way sound sources)
 * Magnetohydrodynamics (MHD)
 * Relativistic Magnetohydrodynamics (RMHD)
-</details>
 
 ### Numerics
 
@@ -157,7 +156,6 @@ They are organized below.
 * RK4-5 operator splitting for Euler-Lagrange modeling
 * Interface sharpening (THINC-like)
 
-
 ### Large-scale and accelerated simulation
 
 * GPU compatible on NVIDIA ([P/V/A/H]100, GH200, etc.) and AMD (MI[1/2/3]00+) GPU and APU hardware
@@ -182,27 +180,20 @@ They are organized below.
 
 ## Citation
 
-If you use MFC, consider citing it as:
-
-<p align="center">
-  <a href="https://doi.org/10.1016/j.cpc.2020.107396">
-    S. H. Bryngelson, K. Schmidmayer, V. Coralic, K. Maeda, J. Meng, T. Colonius (2021) Computer Physics Communications <b>266</b>, 107396
-  </a>
-</p>
+If you use MFC, consider citing it as below.
+Ref. 1 includes all modern MFC features, including GPU acceleration and many new physics features.
+If referencing MFC's (GPU) performance, consider citing ref. 1 and 2, which describe the solver and how it was crafted.
+The original open-source release of MFC is ref. 3, which should be cited for provenance as appropriate.
 
 ```bibtex
-@article{Bryngelson_2021,
-  title   = {{MFC: A}n open-source high-order multi-component, multi-phase, and multi-scale compressible flow solver},
-  author  = {S. H. Bryngelson and K. Schmidmayer and V. Coralic and J. C. Meng and K. Maeda and T. Colonius},
-  journal = {Computer Physics Communications},
-  year    = {2021},
-  volume  = {266},
-  pages   = {107396},
-  doi     = {10.1016/j.cpc.2020.107396}
+@article{Wilfong_2025,
+  author = {Wilfong, Benjamin and {Le Berre}, Henry and Radhakrishnan, Anand and Gupta, Ansh and Vaca-Revelo, Diego and Adam, Dimitrios and Yu, Haocheng and Lee, Hyeoksu and Chreim, Jose Rodolfo and {Carcana Barbosa}, Mirelys and Zhang, Yanjun and Cisneros-Garibay, Esteban and Gnanaskandan, Aswin and {Rodriguez Jr.}, Mauro and Budiardja, Reuben D. and Abbott, Stephen and Colonius, Tim and Bryngelson, Spencer H.},
+  title = {{MFC 5.0: A}n exascale many-physics flow solver},
+  journal = {arXiv preprint arXiv:2503.07953},
+  year = {2025},
+  doi = {10.48550/arXiv.2503.07953}
 }
-```
 
-```bibtex
 @article{Radhakrishnan_2024,
   title   = {Method for portable, scalable, and performant {GPU}-accelerated simulation of multiphase compressible flow},
   author  = {A. Radhakrishnan and H. {Le Berre} and B. Wilfong and J.-S. Spratt and M. {Rodriguez Jr.} and T. Colonius and S. H. Bryngelson},
@@ -212,6 +203,16 @@ If you use MFC, consider citing it as:
   pages   = {109238},
   doi     = {10.1016/j.cpc.2024.109238}
 }
+
+@article{Bryngelson_2021,
+  title   = {{MFC: A}n open-source high-order multi-component, multi-phase, and multi-scale compressible flow solver},
+  author  = {S. H. Bryngelson and K. Schmidmayer and V. Coralic and J. C. Meng and K. Maeda and T. Colonius},
+  journal = {Computer Physics Communications},
+  year    = {2021},
+  volume  = {266},
+  pages   = {107396},
+  doi     = {10.1016/j.cpc.2020.107396}
+}
 ```
 
 ## License
@@ -225,12 +226,7 @@ Federal sponsors have supported MFC development, including the US Department of
 
 MFC computations have used many supercomputing systems. A partial list is below
   * OLCF Frontier and Summit, and testbeds Wombat, Crusher, and Spock (allocation CFD154, PI Bryngelson)
-  * LLNL Tuolumne and Lassen, El Capitan early access system Tioga
-  * PSC Bridges(1/2), NCSA Delta, SDSC Comet and Expanse, Purdue Anvil, TACC Stampede(1-3), and TAMU ACES via ACCESS-CI allocations from Bryngelson, Colonius, Rodriguez, and more.
-  * DOD systems Onyx, Carpenter, Nautilus, and Narwhal via the DOD HPCMP program
+  * LLNL El Capitan, Tuolumne, and Lassen; El Capitan early access system Tioga
+  * NCSA Delta and DeltaAI, PSC Bridges(1/2), SDSC Comet and Expanse, Purdue Anvil, TACC Stampede(1-3), and TAMU ACES via ACCESS-CI allocations from Bryngelson, Colonius, Rodriguez, and more.
+  * DOD systems Blueback, Onyx, Carpenter, Nautilus, and Narwhal via the DOD HPCMP program
   * Sandia National Labs systems Doom and Attaway and testbed systems Weaver and Vortex
-
-
-## Contributors
-
-[![Contributors](https://contributors-img.web.app/image?repo=mflowcode/mfc)](https://github.com/mflowcode/mfc/graphs/contributors)

docs/documentation/readme.md

@@ -9,6 +9,7 @@
 - [Running MFC](running.md)
 - [Flow Visualization](visualization.md)
 - [Performance](expectedPerformance.md)
+- [GPU Parallelization](gpuParallelization.md)
 - [MFC's Authors](authors.md)
 - [References](references.md)
 

src/common/include/macros.fpp

@@ -96,9 +96,6 @@
     end if
 #:enddef
 
-#define t_vec3   real(wp), dimension(1:3)
-#define t_mat4x4 real(wp), dimension(1:4,1:4)
-
 #:def ASSERT(predicate, message = None)
     if (.not. (${predicate}$)) then
         call s_mpi_abort("${_FILE_.split('/')[-1]}$:${_LINE_}$: "// &

src/common/m_derived_types.fpp

@@ -141,13 +141,13 @@ module m_derived_types
         character(LEN=pathlen_max) :: filepath !<
         !! Path the STL file relative to case_dir.
 
-        t_vec3 :: translate !<
+        real(wp), dimension(1:3) :: translate !<
         !! Translation of the STL object.
 
-        t_vec3 :: scale !<
+        real(wp), dimension(1:3) :: scale !<
         !! Scale factor for the STL object.
 
-        t_vec3 :: rotate !<
+        real(wp), dimension(1:3) :: rotate !<
         !! Angle to rotate the STL object along each cartesian coordinate axis,
         !! in radians.
 
@@ -160,17 +160,17 @@ module m_derived_types
 
     type :: t_triangle
         real(wp), dimension(1:3, 1:3) :: v ! Vertices of the triangle
-        t_vec3 :: n ! Normal vector
+        real(wp), dimension(1:3) :: n ! Normal vector
     end type t_triangle
 
     type :: t_ray
-        t_vec3 :: o ! Origin
-        t_vec3 :: d ! Direction
+        real(wp), dimension(1:3) :: o ! Origin
+        real(wp), dimension(1:3) :: d ! Direction
     end type t_ray
 
     type :: t_bbox
-        t_vec3 :: min ! Minimum coordinates
-        t_vec3 :: max ! Maximum coordinates
+        real(wp), dimension(1:3) :: min ! Minimum coordinates
+        real(wp), dimension(1:3) :: max ! Maximum coordinates
     end type t_bbox
 
     type :: t_model
@@ -265,13 +265,13 @@ module m_derived_types
         character(LEN=pathlen_max) :: model_filepath !<
         !! Path the STL file relative to case_dir.
 
-        t_vec3 :: model_translate !<
+        real(wp), dimension(1:3) :: model_translate !<
         !! Translation of the STL object.
 
-        t_vec3 :: model_scale !<
+        real(wp), dimension(1:3) :: model_scale !<
         !! Scale factor for the STL object.
 
-        t_vec3 :: model_rotate !<
+        real(wp), dimension(1:3) :: model_rotate !<
         !! Angle to rotate the STL object along each cartesian coordinate axis,
         !! in radians.
 
@@ -303,13 +303,13 @@ module m_derived_types
         character(LEN=pathlen_max) :: model_filepath !<
         !! Path the STL file relative to case_dir.
 
-        t_vec3 :: model_translate !<
+        real(wp), dimension(1:3) :: model_translate !<
         !! Translation of the STL object.
 
-        t_vec3 :: model_scale !<
+        real(wp), dimension(1:3) :: model_scale !<
         !! Scale factor for the STL object.
 
-        t_vec3 :: model_rotate !<
+        real(wp), dimension(1:3) :: model_rotate !<
         !! Angle to rotate the STL object along each cartesian coordinate axis,
         !! in radians.
 
@@ -443,4 +443,10 @@ module m_derived_types
 
     end type bubbles_lagrange_parameters
 
+    !> Max and min number of cells in a direction of each combination of x-,y-, and z-
+    type cell_num_bounds
+        integer :: mn_max, np_max, mp_max, mnp_max
+        integer :: mn_min, np_min, mp_min, mnp_min
+    end type cell_num_bounds
+
 end module m_derived_types

src/common/m_helper.fpp

@@ -319,8 +319,8 @@ contains
     pure function f_create_transform_matrix(p, center) result(out_matrix)
 
         type(ic_model_parameters), intent(in) :: p
-        t_vec3, optional, intent(in) :: center
-        t_mat4x4 :: sc, rz, rx, ry, tr, t_back, t_to_origin, out_matrix
+        real(wp), dimension(1:3), optional, intent(in) :: center
+        real(wp), dimension(1:4, 1:4) :: sc, rz, rx, ry, tr, t_back, t_to_origin, out_matrix
 
         sc = transpose(reshape([ &
                                p%scale(1), 0._wp, 0._wp, 0._wp, &
@@ -379,8 +379,8 @@ contains
     !! @param matrix Transformation matrix.
     pure subroutine s_transform_vec(vec, matrix)
 
-        t_vec3, intent(inout) :: vec
-        t_mat4x4, intent(in) :: matrix
+        real(wp), dimension(1:3), intent(inout) :: vec
+        real(wp), dimension(1:4, 1:4), intent(in) :: matrix
 
         real(wp), dimension(1:4) :: tmp
 
@@ -395,7 +395,7 @@ contains
     pure subroutine s_transform_triangle(triangle, matrix, matrix_n)
 
         type(t_triangle), intent(inout) :: triangle
-        t_mat4x4, intent(in) :: matrix, matrix_n
+        real(wp), dimension(1:4, 1:4), intent(in) :: matrix, matrix_n
 
         integer :: i
 
@@ -413,7 +413,7 @@ contains
     pure subroutine s_transform_model(model, matrix, matrix_n)
 
         type(t_model), intent(inout) :: model
-        t_mat4x4, intent(in) :: matrix, matrix_n
+        real(wp), dimension(1:4, 1:4), intent(in) :: matrix, matrix_n
 
         integer :: i
 

src/common/m_helper_basic.fpp

@@ -16,7 +16,8 @@ module m_helper_basic
               f_is_default, &
               f_all_default, &
               f_is_integer, &
-              s_configure_coordinate_bounds
+              s_configure_coordinate_bounds, &
+              s_update_cell_bounds
 
 contains
 
@@ -146,4 +147,24 @@ contains
 
     end subroutine s_configure_coordinate_bounds
 
+    !> Updates the min and max number of cells in each set of axes
+    !! @param bounds Min ans max values to update
+    !! @param m Number of cells in x-axis
+    !! @param n Number of cells in y-axis
+    !! @param p Number of cells in z-axis
+    pure elemental subroutine s_update_cell_bounds(bounds, m, n, p)
+        type(cell_num_bounds), intent(out) :: bounds
+        integer, intent(in) :: m, n, p
+
+        bounds%mn_max = max(m, n)
+        bounds%np_max = max(n, p)
+        bounds%mp_max = max(m, p)
+        bounds%mnp_max = max(m, n, p)
+        bounds%mn_min = min(m, n)
+        bounds%np_min = min(n, p)
+        bounds%mp_min = min(m, p)
+        bounds%mnp_min = min(m, n, p)
+
+    end subroutine s_update_cell_bounds
+
 end module m_helper_basic

src/common/m_mpi_common.fpp

@@ -67,10 +67,10 @@ contains
                                          & (m + 2*buff_size + 1)* &
                                          & (n + 2*buff_size + 1)* &
                                          & (p + 2*buff_size + 1)/ &
-                                         & (min(m, n, p) + 2*buff_size + 1))
+                                         & (cells_bounds%mnp_min + 2*buff_size + 1))
             else
                 halo_size = -1 + buff_size*(v_size)* &
-                                         & (max(m, n) + 2*buff_size + 1)
+                                         & (cells_bounds%mn_max + 2*buff_size + 1)
             end if
         else
             halo_size = -1 + buff_size*(v_size)
@@ -1446,6 +1446,8 @@ contains
             end if
         end do
 
+        call s_update_cell_bounds(cells_bounds, m, n, p)
+
         ! Boundary condition at the beginning
         if (proc_coords(1) > 0 .or. (bc_x%beg == BC_PERIODIC .and. num_procs_x > 1)) then
             proc_coords(1) = proc_coords(1) - 1

src/common/m_variables_conversion.fpp

@@ -637,7 +637,7 @@ contains
 #ifdef MFC_SIMULATION
 
         if (viscous) then
-            @:ALLOCATE(Res(1:2, 1:maxval(Re_size)))
+            @:ALLOCATE(Res(1:2, 1:Re_size_max))
             do i = 1, 2
                 do j = 1, Re_size(i)
                     Res(i, j) = fluid_pp(Re_idx(i, j))%Re(i)

src/post_process/m_global_parameters.fpp

@@ -38,6 +38,9 @@ module m_global_parameters
     integer :: p
     !> @}
 
+    !> @name Max and min number of cells in a direction of each combination of x-,y-, and z-
+    type(cell_num_bounds) :: cells_bounds
+
     integer(8) :: nGlobal ! Total number of cells in global domain
 
     !> @name Cylindrical coordinates (either axisymmetric or full 3D)
@@ -336,6 +339,8 @@ contains
 
         ! Computational domain parameters
         m = dflt_int; n = 0; p = 0
+        call s_update_cell_bounds(cells_bounds, m, n, p)
+
         m_root = dflt_int
         cyl_coord = .false.
 

src/post_process/m_start_up.f90

@@ -111,6 +111,9 @@ impure subroutine s_read_input_file
             end if
 
             close (1)
+
+            call s_update_cell_bounds(cells_bounds, m, n, p)
+
             ! Store m,n,p into global m,n,p
             m_glb = m
             n_glb = n