Extend initial conditions from 1D to 2D and 2D to 3D (#844)

DimAdam-01 · sbryngelson · web-flow · commit 7bb8ed31e564 · 2025-06-21T12:19:29.000-04:00
Co-authored-by: Spencer Bryngelson &lt;sbryngelson@gmail.com&gt;
diff --git a/docs/documentation/case.md b/docs/documentation/case.md
@@ -236,6 +236,16 @@ and use `patch_icpp(i)%%geometry = 7` and `patch_icpp(i)%%hcid = 200` in the inp
 Additional variables can be declared in `Hardcoded1[2,3]DVariables` and used in `hardcoded1[2,3]D`.
 As a convention, any hard coded patches that are part of the MFC master branch should be identified as 1[2,3]xx where the first digit indicates the number of dimensions.
 
+The code provides three pre-built patches for dimensional extrusion of initial conditions:
+
+- `case(170)`: Load 1D profile from data files
+- `case(270)`: Extrude 1D data to 2D domain
+- `case(370)`: Extrude 2D data to 3D domain
+
+Setup: Only requires specifying `init_dir` and filename pattern via `zeros_default`. Grid dimensions are automatically detected from the data files.
+Implementation: All variables and file handling are managed in `src/pre_process/include/ExtrusionHardcodedIC.fpp` with no manual grid configuration needed.
+Usage: Ideal for initializing simulations from lower-dimensional solutions, enabling users to add perturbations or modifications to the base extruded fields for flow instability studies.
+
 #### Parameter Descriptions
 
 - `num_patches` defines the total number of patches defined in the domain.
diff --git a/src/pre_process/include/1dHardcodedIC.fpp b/src/pre_process/include/1dHardcodedIC.fpp
@@ -1,13 +1,14 @@
 #:def Hardcoded1DVariables()
     ! Place any declaration of intermediate variables here
-
 #:enddef
 
 #:def Hardcoded1D()
 
     select case (patch_icpp(patch_id)%hcid)
-    case (100)
-        ! Put your variable assignments here
+    case (170)
+        ! This hardcoded case can be used to start a simulation with initial conditions given from a known 1D profile (e.g. Cantera, SDtoolbox)
+        @: HardcodedReadValues()
+
     case default
         call s_int_to_str(patch_id, iStr)
         call s_mpi_abort("Invalid hcid specified for patch "//trim(iStr))
diff --git a/src/pre_process/include/2dHardcodedIC.fpp b/src/pre_process/include/2dHardcodedIC.fpp
@@ -1,12 +1,11 @@
 #:def Hardcoded2DVariables()
-
+    ! Place any declaration of intermediate variables here
     real(wp) :: eps
     real(wp) :: r, rmax, gam, umax, p0
     real(wp) :: rhoH, rhoL, pRef, pInt, h, lam, wl, amp, intH, intL, alph
     real(wp) :: factor
-
+    
     eps = 1.e-9_wp
-
 #:enddef
 
 #:def Hardcoded2D()
@@ -158,6 +157,10 @@
             q_prim_vf(E_idx)%sf(i, j, 0) = 3.e-5_wp
         end if
 
+    case (270)
+        ! This hardcoded case extrudes a 1D profile to initialize a 2D simulation domain
+        @: HardcodedReadValues()
+
     case default
         if (proc_rank == 0) then
             call s_int_to_str(patch_id, iStr)
diff --git a/src/pre_process/include/3dHardcodedIC.fpp b/src/pre_process/include/3dHardcodedIC.fpp
@@ -1,6 +1,5 @@
 #:def Hardcoded3DVariables()
     ! Place any declaration of intermediate variables here
-
     real(wp) :: rhoH, rhoL, pRef, pInt, h, lam, wl, amp, intH, alph
 
     real(wp) :: eps
@@ -56,7 +55,10 @@
             q_prim_vf(advxe)%sf(i, j, k) = patch_icpp(1)%alpha(2)
         end if
 
-        ! Put your variable assignments here
+    case (370)
+        ! This hardcoded case extrudes a 2D profile to initialize a 3D simulation domain
+        @: HardcodedReadValues()
+
     case default
         call s_int_to_str(patch_id, iStr)
         call s_mpi_abort("Invalid hcid specified for patch "//trim(iStr))
diff --git a/src/pre_process/include/ExtrusionHardcodedIC.fpp b/src/pre_process/include/ExtrusionHardcodedIC.fpp
@@ -0,0 +1,209 @@
+!> @brief Allocate memory and read initial condition data for IC extrusion.
+!>
+!> @details
+!>   This macro handles the complete initialization process for IC extrusion by:
+!>
+!>   **Memory Allocation:**
+!>     - stored_values(xRows, yRows, sys_size) - stores primitive variable data from files
+!>     - x_coords(nrows) - stores x-coordinates from input files
+!>     - y_coords(nrows) - stores y-coordinates from input files (3D case only)
+!>
+!>   **File Reading Operations:**
+!>     - Reads primitive variable data from multiple files with pattern:
+!>       `prim.<file_number>.00.<timestep>.dat` where timestep uses `zeros_default` padding
+!>     - Files are read from directory specified by `init_dir` parameter
+!>     - Supports 1D, 2D, and 3D computational domains
+!>
+!>   **Grid Structure Detection:**
+!>     - 1D/2D: Counts lines in first file to determine xRows
+!>     - 3D: Analyzes coordinate patterns to determine xRows and yRows structure
+!>
+!>   **MPI Domain Mapping:**
+!>     - Calculates global_offset_x and global_offset_y for MPI subdomain positioning
+!>     - Maps file coordinates to local computational grid coordinates
+!>
+!>   **Data Assignment:**
+!>     - Populates q_prim_vf primitive variable arrays with file data
+!>     - Handles momentum component indexing with special treatment for momxe
+!>     - Sets momxe component to zero for 2D/3D cases
+!>
+!>   **State Management:**
+!>     - Uses files_loaded flag to prevent redundant file operations
+!>     - Preserves data across multiple macro calls within same simulation
+!>
+!>   @note File pattern uses `zeros_default` parameter (default: "000000") for timestep padding
+!>   @note Directory path is hardcoded in `init_dir` parameter - modify as needed
+!>   @warning Aborts execution if file reading errors occur.
+
+#:def HardcodedDimensionsExtrusion()
+    integer :: xRows, yRows, nRows, iix, iiy, max_files
+    integer :: f, iter, ios, ios2, unit, unit2, idx, idy, index_x, index_y, jump, line_count, ycount
+    real(wp) :: x_len, x_step, y_len, y_step
+    real(wp) :: dummy_x, dummy_y, dummy_z, x0, y0
+    integer :: global_offset_x, global_offset_y           ! MPI subdomain offset
+    real(wp) :: delta_x, delta_y
+    character(len=100), dimension(sys_size) :: fileNames ! Arrays to store all data from files
+    character(len=200) :: errmsg
+    real(wp), allocatable :: stored_values(:, :, :)
+    real(wp), allocatable :: x_coords(:), y_coords(:)
+    logical :: files_loaded = .false.
+    real(wp) :: domain_xstart, domain_xend, domain_ystart, domain_yend
+    character(len=*), parameter :: init_dir = "/home/MFC/FilesDirectory" ! For example /home/MFC/examples/1D_Shock/D/
+    character(len=20) :: file_num_str     ! For storing the file number as a string
+    character(len=20) :: zeros_part       ! For the trailing zeros part
+    character(len=6), parameter :: zeros_default = "000000"  ! Default zeros (can be changed)
+#:enddef
+
+#:def HardcodedReadValues()
+
+    if (.not. files_loaded) then
+        max_files = merge(sys_size, sys_size - 1, num_dims == 1)
+        do f = 1, max_files
+            write (file_num_str, '(I0)') f
+            fileNames(f) = trim(init_dir)//"prim."//trim(file_num_str)//".00."//zeros_default//".dat"
+        end do
+
+        ! Common file reading setup
+        open (newunit=unit2, file=trim(fileNames(1)), status='old', action='read', iostat=ios2)
+        if (ios2 /= 0) call s_mpi_abort("Error opening file: "//trim(fileNames(1)))
+
+        select case (num_dims)
+        case (1, 2)  ! 1D and 2D cases are similar
+            ! Count lines
+            line_count = 0
+            do
+                read (unit2, *, iostat=ios2) dummy_x, dummy_y
+                if (ios2 /= 0) exit
+                line_count = line_count + 1
+            end do
+            close (unit2)
+
+            xRows = line_count
+            yRows = 1
+            index_x = 0
+            if (num_dims == 2) index_x = i
+            @:ALLOCATE (x_coords(xRows), stored_values(xRows, 1, sys_size))
+
+            ! Read data from all files
+            do f = 1, max_files
+                open (newunit=unit, file=trim(fileNames(f)), status='old', action='read', iostat=ios)
+                if (ios /= 0) call s_mpi_abort("Error opening file: "//trim(fileNames(f)))
+
+                do iter = 1, xRows
+                    read (unit, *, iostat=ios) x_coords(iter), stored_values(iter, 1, f)
+                    if (ios /= 0) call s_mpi_abort("Error reading file: "//trim(fileNames(f)))
+                end do
+                close (unit)
+            end do
+
+            ! Calculate offsets
+            domain_xstart = x_coords(1)
+            x_step = x_cc(1) - x_cc(0)
+            delta_x = merge(x_cc(0) - domain_xstart + x_step/2.0, &
+                            x_cc(index_x) - domain_xstart + x_step/2.0, num_dims == 1)
+            global_offset_x = nint(abs(delta_x)/x_step)
+
+        case (3)  ! 3D case - determine grid structure
+            ! Find yRows by counting rows with same x
+            read (unit2, *, iostat=ios2) x0, y0, dummy_z
+            if (ios2 /= 0) call s_mpi_abort("Error reading first line")
+
+            yRows = 1
+            do
+                read (unit2, *, iostat=ios2) dummy_x, dummy_y, dummy_z
+                if (ios2 /= 0) exit
+                if (dummy_x == x0 .and. dummy_y /= y0) then
+                    yRows = yRows + 1
+                else
+                    exit
+                end if
+            end do
+            close (unit2)
+
+            ! Count total rows
+            open (newunit=unit2, file=trim(fileNames(1)), status='old', action='read', iostat=ios2)
+            nrows = 0
+            do
+                read (unit2, *, iostat=ios2) dummy_x, dummy_y, dummy_z
+                if (ios2 /= 0) exit
+                nrows = nrows + 1
+            end do
+            close (unit2)
+
+            xRows = nrows/yRows
+            @:ALLOCATE (x_coords(nrows), y_coords(nrows), stored_values(xRows, yRows, sys_size))
+            index_x = i
+            index_y = j
+
+            ! Read all files
+            do f = 1, max_files
+                open (newunit=unit, file=trim(fileNames(f)), status='old', action='read', iostat=ios)
+                if (ios /= 0) then
+                    if (f == 1) call s_mpi_abort("Error opening file: "//trim(fileNames(f)))
+                    cycle
+                end if
+
+                iter = 0
+                do iix = 1, xRows
+                    do iiy = 1, yRows
+                        iter = iter + 1
+                        if (f == 1) then
+                            read (unit, *, iostat=ios) x_coords(iter), y_coords(iter), stored_values(iix, iiy, f)
+                        else
+                            read (unit, *, iostat=ios) dummy_x, dummy_y, stored_values(iix, iiy, f)
+                        end if
+                        if (ios /= 0) call s_mpi_abort("Error reading data")
+                    end do
+                end do
+                close (unit)
+            end do
+
+            ! Calculate offsets
+            x_step = x_cc(1) - x_cc(0)
+            y_step = y_cc(1) - y_cc(0)
+            delta_x = x_cc(index_x) - x_coords(1) + x_step/2.0_wp
+            delta_y = y_cc(index_y) - y_coords(1) + y_step/2.0_wp
+            global_offset_x = nint(abs(delta_x)/x_step)
+            global_offset_y = nint(abs(delta_y)/y_step)
+        end select
+
+        files_loaded = .true.
+    end if
+
+    ! Data assignment
+    select case (num_dims)
+    case (1)
+        idx = i + 1 + global_offset_x
+        do f = 1, sys_size
+            q_prim_vf(f)%sf(i, 0, 0) = stored_values(idx, 1, f)
+        end do
+
+    case (2)
+        idx = i + 1 + global_offset_x - index_x
+        do f = 1, sys_size - 1
+            jump = merge(1, 0, f >= momxe)
+            q_prim_vf(f + jump)%sf(i, j, 0) = stored_values(idx, 1, f)
+        end do
+        q_prim_vf(momxe)%sf(i, j, 0) = 0.0_wp
+
+    case (3)
+        idx = i + 1 + global_offset_x - index_x
+        idy = j + 1 + global_offset_y - index_y
+        do f = 1, sys_size - 1
+            jump = merge(1, 0, f >= momxe)
+            q_prim_vf(f + jump)%sf(i, j, k) = stored_values(idx, idy, f)
+        end do
+        q_prim_vf(momxe)%sf(i, j, k) = 0.0_wp
+    end select
+#:enddef
+
+#:def HardcodedDellacation()
+    if (allocated(stored_values)) then
+        @:DEALLOCATE (stored_values)
+        @:DEALLOCATE (x_coords)
+    end if
+
+    if (allocated(y_coords)) then
+        @:DEALLOCATE (y_coords)
+    end if
+#:enddef
diff --git a/src/pre_process/m_patches.fpp b/src/pre_process/m_patches.fpp
@@ -3,6 +3,7 @@
 !! @brief Contains module m_patches
 
 #:include 'case.fpp'
+#:include 'ExtrusionHardcodedIC.fpp'
 #:include '1dHardcodedIC.fpp'
 #:include '2dHardcodedIC.fpp'
 #:include '3dHardcodedIC.fpp'
@@ -1296,15 +1297,18 @@ contains
         type(scalar_field), dimension(1:sys_size), intent(inout) :: q_prim_vf
 
         ! Generic loop iterators
-        integer :: i
+        integer :: i, j, k
+
         ! Placeholders for the cell boundary values
         real(wp) :: pi_inf, gamma, lit_gamma
-
+        @:HardcodedDimensionsExtrusion()
         @:Hardcoded1DVariables()
 
         pi_inf = fluid_pp(1)%pi_inf
         gamma = fluid_pp(1)%gamma
         lit_gamma = (1._wp + gamma)/gamma
+        j = 0.0_wp
+        k = 0.0_wp
 
         ! Transferring the patch's centroid and length information
         x_centroid = patch_icpp(patch_id)%x_centroid
@@ -1341,6 +1345,8 @@ contains
             end if
         end do
 
+        @:HardcodedDellacation()
+
     end subroutine s_1D_analytical
 
         !! @param patch_id is the patch identifier
@@ -1409,17 +1415,19 @@ contains
         integer, dimension(0:m, 0:n, 0:p), intent(inout) :: patch_id_fp
         type(scalar_field), dimension(1:sys_size), intent(inout) :: q_prim_vf
 
-        integer :: i, j !< generic loop iterators
+        integer :: i, j, k !< generic loop iterators
 
         real(wp) :: pi_inf, gamma, lit_gamma !< equation of state parameters
         real(wp) :: l, U0 !< Taylor Green Vortex parameters
-
+        @:HardcodedDimensionsExtrusion()
         @:Hardcoded2DVariables()
 
         pi_inf = fluid_pp(1)%pi_inf
         gamma = fluid_pp(1)%gamma
         lit_gamma = (1._wp + gamma)/gamma
 
+        k = 0.0_wp
+
         ! Transferring the patch's centroid and length information
         x_centroid = patch_icpp(patch_id)%x_centroid
         y_centroid = patch_icpp(patch_id)%y_centroid
@@ -1465,6 +1473,8 @@ contains
             end do
         end do
 
+        @:HardcodedDellacation()
+
     end subroutine s_2D_analytical
 
     !> This patch assigns the primitive variables as analytical
@@ -1480,7 +1490,7 @@ contains
 
         integer :: i, j, k !< generic loop iterators
         real(wp) :: pi_inf, gamma, lit_gamma !< equation of state parameters
-
+        @:HardcodedDimensionsExtrusion()
         @:Hardcoded3DVariables()
 
         pi_inf = fluid_pp(1)%pi_inf
@@ -1549,6 +1559,7 @@ contains
             end do
         end do
 
+        @:HardcodedDellacation()
     end subroutine s_3D_analytical
 
     !> This patch generates the shape of the spherical harmonics
diff --git a/src/simulation/m_data_output.fpp b/src/simulation/m_data_output.fpp
@@ -663,8 +663,6 @@ contains
                             if (((i >= cont_idx%beg) .and. (i <= cont_idx%end)) &
                                 .or. &
                                 ((i >= adv_idx%beg) .and. (i <= adv_idx%end)) &
-                                .or. &
-                                ((i >= chemxb) .and. (i <= chemxe)) &
                                 ) then
                                 write (2, FMT) x_cb(j), y_cb(k), q_cons_vf(i)%sf(j, k, 0)
                             else
