Extend initial conditions from 1D to 2D and 2D to 3D

docs/documentation/case.md

@@ -236,6 +236,11 @@ 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.
 
+Also Several custom hard-coded patches are already defined to support common workflows. These include:
+
+By convention, patch ID `case(170)` load a 1D profile, ID `case(270)` extrude that 1D data into 2D, and ID `case(370)` extrude 2D data into 3D. 
+You only need to supply `init_dir` (and the filename‐pattern via `zeros_default`). All related variables (including init_dir, zeros_default, and the “files_loaded” flag) 
+live in `src/pre_process/include/ExtrusionHardcodedIC.fpp`, and no manual grid‐size settings are required.
 #### Parameter Descriptions
 
 - `num_patches` defines the total number of patches defined in the domain.

src/pre_process/include/1dHardcodedIC.fpp

@@ -1,13 +1,75 @@
 #: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)
+
+        ! Generate file names in a loop
+        do f = 1, sys_size
+            ! Convert file number to string with proper formatting
+            if (f < 10) then
+                write (file_num_str, '(I1)') f  ! Single digit
+            else
+                write (file_num_str, '(I2)') f  ! Double digit
+            end if
+            fileNames(f) = trim(init_dir)//"prim."//trim(file_num_str)//".00."//zeros_default//".dat"
+            ! Create the filename with the pattern "prim.X.00.000000.dat"
+        end do
+
+        if (.not. files_loaded) then
+            !Reading the first file to calculate the number of grid points in the x direction
+            line_count = 0
+            open (newunit=unit2, file=trim(fileNames(1)), status='old', action='read', iostat=ios2)
+            if (ios2 /= 0) then
+                write (errmsg, '(A,A)') "Error opening file: ", trim(fileNames(1))
+                call s_mpi_abort(trim(errmsg))
+            end if
+            do
+                read (unit2, *, iostat=ios2) dummy_x, dummy_y
+                if (ios2 /= 0) exit  ! Exit since the file has been read
+                line_count = line_count + 1
+            end do
+            close (unit2)
+
+            xRows = line_count
+            allocate (x_coords(xRows))
+            allocate (stored_values(xRows, 1, sys_size))
+            do f = 1, sys_size
+                ! Open the file for reading
+                open (newunit=unit, file=trim(fileNames(f)), status='old', action='read', iostat=ios)
+                if (ios /= 0) then
+                    write (errmsg, '(A,A)') "Error opening file: ", trim(fileNames(f))
+                    call s_mpi_abort(trim(errmsg))
+                end if
+                ! Read all rows at once into memory
+                do iter = 1, xRows
+                    read (unit, *, iostat=ios) x_coords(iter), stored_values(iter, 1, f)
+                    if (ios /= 0) then
+                        write (errmsg, '(A,A,A,I0,A)') ' Error reading "', trim(fileNames(f)), &
+                            '" at index (', iter, ')'
+                        call s_mpi_abort(trim(errmsg)) ! Exit loop on error
+                    end if
+                end do
+                close (unit)
+            end do
+            ! Calculate domain information for mapping
+            domain_xstart = x_coords(1)
+            domain_xend = x_coords(xRows)
+            x_step = x_cc(1) - x_cc(0)
+
+            delta_x = x_cc(0) - domain_xstart - x_step/2.0
+            global_offset_x = nint(abs(delta_x)/x_step)
+            files_loaded = .true.
+        end if
+        ! Simple mapping - find the closest index
+        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 default
         call s_int_to_str(patch_id, iStr)
         call s_mpi_abort("Invalid hcid specified for patch "//trim(iStr))

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 = 1e-9_wp
-
 #:enddef
 
 #:def Hardcoded2D()
@@ -158,6 +157,77 @@
             q_prim_vf(E_idx)%sf(i, j, 0) = 3.e-5_wp
         end if
 
+    case (270)
+
+        do f = 1, sys_size - 1
+            ! Convert file number to string with proper formatting
+            if (f < 10) then
+                write (file_num_str, '(I1)') f  ! Single digit
+            else
+                write (file_num_str, '(I2)') f  ! Double digit
+                ! For more than 99 files, you might need to adjust this format
+            end if
+            fileNames(f) = trim(init_dir)//"prim."//trim(file_num_str)//".00."//zeros_default//".dat"
+        end do
+
+        if (.not. files_loaded) then
+            ! Calculating the number of grid points in the x direction in the file.
+            line_count = 0
+            open (newunit=unit2, file=trim(fileNames(1)), status='old', action='read', iostat=ios2)
+            if (ios2 /= 0) then
+                write (errmsg, '(A,A)') "Error opening file: ", trim(fileNames(1))
+                call s_mpi_abort(trim(errmsg))
+            end if
+            do
+                read (unit2, *, iostat=ios2) dummy_x, dummy_y
+                if (ios2 /= 0) exit ! Exit since files has been read
+                line_count = line_count + 1
+            end do
+            close (unit2)
+
+            index_x = i
+            xRows = line_count
+            allocate (x_coords(xRows))
+            allocate (stored_values(xRows, 1, sys_size))
+            do f = 1, sys_size - 1
+                ! Open the file for reading
+                open (newunit=unit, file=trim(fileNames(f)), status='old', action='read', iostat=ios)
+                if (ios /= 0 .and. proc_rank == 0) then
+                    write (errmsg, '(A,A)') "Error opening file: ", trim(fileNames(f))
+                    call s_mpi_abort(trim(errmsg))
+                end if
+                ! Read all rows at once into memory
+                do iter = 1, xRows
+                    read (unit, *, iostat=ios) x_coords(iter), stored_values(iter, 1, f)
+                    if (ios /= 0) then
+                        write (errmsg, '(A,A,A,I0,A)') 'Error reading "', trim(fileNames(f)), &
+                            '" at index (', iter, ')'
+                        call s_mpi_abort(trim(errmsg)) ! Exit loop on error
+                    end if
+                end do
+                close (unit)
+            end do
+            ! Calculate domain information for mapping
+            domain_xstart = x_coords(1)
+            domain_xend = x_coords(xRows)
+            x_step = x_cc(1) - x_cc(0)
+            delta_x = x_cc(index_x) - domain_xstart + x_step/2
+            global_offset_x = nint(abs(delta_x)/x_step)
+            files_loaded = .true.
+        end if
+        ! Calculate the index in the file data array corresponding to x_cc(i)
+        idx = i + 1 + global_offset_x - index_x
+        do f = 1, sys_size - 1
+            if (f >= momxe) then
+                jump = 1
+            else
+                jump = 0
+            end if
+            q_prim_vf(f + jump)%sf(i, j, 0) = stored_values(idx, 1, f)
+        end do
+        ! Set element the velocity paraller to y explicitly to zero
+        q_prim_vf(momxe)%sf(i, j, 0) = 0.0_wp
+
     case default
         if (proc_rank == 0) then
             call s_int_to_str(patch_id, iStr)

src/pre_process/include/3dHardcodedIC.fpp

@@ -1,10 +1,8 @@
 #:def Hardcoded3DVariables()
     ! Place any declaration of intermediate variables here
-
     real(wp) :: rhoH, rhoL, pRef, pInt, h, lam, wl, amp, intH, alph
 
     real(wp) :: eps
-
     eps = 1e-9_wp
 #:enddef
 
@@ -56,6 +54,150 @@
             q_prim_vf(advxe)%sf(i, j, k) = patch_icpp(1)%alpha(2)
         end if
 
+    case (370)
+        ! Generate file names dynamically in a loop
+        do f = 1, sys_size - 1
+            ! Convert file number to string with proper formatting
+            if (f < 10) then
+                write (file_num_str, '(I1)') f  ! Single digit
+            else
+                write (file_num_str, '(I2)') f  ! Double digit
+                ! For more than 99 files, you might need to adjust this format
+            end if
+            fileNames(f) = trim(init_dir)//"prim."//trim(file_num_str)//".00."//zeros_default//".dat"
+        end do
+
+        if (.not. files_loaded) then
+            ! 1) Discover yRows by reading fileNames(1)  until x changes from the first line’s x.
+            open (newunit=unit2, file=trim(fileNames(1)), status='old', action='read', iostat=ios2)
+            if (ios2 /= 0) then
+                write (errmsg, '(A,A)') "Error opening file for yRows detection: ", trim(fileNames(1))
+                call s_mpi_abort(trim(errmsg))
+            end if
+
+            ! Read the very first line to get (x₀, y₀)
+            read (unit2, *, iostat=ios2) x0, y0, dummy_z
+            if (ios2 /= 0) then
+                write (errmsg, '(A,A)') "Error reading first line of: ", trim(fileNames(1))
+                call s_mpi_abort(trim(errmsg))
+            end if
+
+            ycount = 1
+            do
+                read (unit2, *, iostat=ios2) dummy_x, dummy_y, dummy_z
+                if (ios2 /= 0) exit         ! End of File : stop counting
+                if (dummy_x == x0 .and. dummy_y /= y0) then
+                    ! As soon as x or y changes from the first line’s (x0,y0),
+                    ! we know we have counted all the (x0, y1) rows.
+                    ycount = ycount + 1
+                else
+                    exit
+                end if
+            end do
+            yRows = ycount
+            close (unit2)
+            ! 2) Count total rows (nrows) to get xRows
+            open (newunit=unit2, file=trim(fileNames(1)), status='old', action='read', iostat=ios2)
+            if (ios2 /= 0) then
+                write (errmsg, '(A,A)') "Error re‐opening file to count rows: ", trim(fileNames(1))
+                call s_mpi_abort(trim(errmsg))
+            end if
+
+            nrows = 0
+            do
+                read (unit2, *, iostat=ios2) dummy_x, dummy_y, dummy_z
+                if (ios2 /= 0) exit
+                nrows = nrows + 1
+            end do
+            close (unit2)
+
+            if (mod(nrows, yRows) /= 0) then
+                write (errmsg, '(A,A,I0,A,I0)') &
+                    "File ’", trim(fileNames(1)), &
+                    "’ has total lines=", nrows, &
+                    " which is not a multiple of yRows=", yRows
+                call s_mpi_abort(trim(errmsg))
+            end if
+            xRows = nrows/yRows
+            allocate (x_coords(nRows))
+            allocate (y_coords(nRows))
+            allocate (stored_values(xRows, yRows, sys_size))
+            index_x = i
+            index_y = j
+
+            open (newunit=unit, file=trim(fileNames(1)), status='old', action='read', iostat=ios)
+            if (ios /= 0) then
+                write (errmsg, '(A,A)') "Error opening file: ", trim(fileNames(f))
+                call s_mpi_abort(trim(errmsg))
+            end if
+
+            iter = 0
+            do iix = 1, xRows
+                do iiy = 1, yRows
+                    iter = iter + 1
+                    read (unit, *, iostat=ios) x_coords(iter), y_coords(iter), stored_values(iix, iiy, 1)
+                    if (ios /= 0) then
+                        write (errmsg, '(A,A,A,I0,A)') 'Error reading "', trim(fileNames(1)), &
+                            '" at indices (', iix, ',', iiy, ')'
+                        call s_mpi_abort(trim(errmsg))
+                    end if
+                end do
+            end do
+            close (unit)
+
+            !Now read only the values from remaining files (skip x,y columns)
+            do f = 2, sys_size - 1
+                open (newunit=unit, file=trim(fileNames(f)), status='old', action='read', iostat=ios)
+                if (ios /= 0) then
+                    print *, "Error opening file: ", trim(fileNames(f))
+                    cycle  ! Skip this file on error
+                end if
+                do iix = 1, xRows
+                    do iiy = 1, yRows
+                        ! Read and discard x,y, then read the value
+                        read (unit, *, iostat=ios) dummy_x, dummy_y, stored_values(iix, iiy, f)
+                        if (ios /= 0) then
+                            write (errmsg, '(A,A,I0,A,I0,A)') 'Error reading "', trim(fileNames(f)), &
+                                '" at indices (', iix, ',', iiy, ')'
+                            call s_mpi_abort(trim(errmsg))
+                        end if
+                    end do
+                end do
+            end do
+            close (unit)
+
+            domain_xstart = x_coords(1)  ! First x value
+            domain_ystart = y_coords(1)  ! First y value
+            ! Calculate actual domain end values
+            domain_xend = x_coords(nRows - yRows + 1)  ! Last x value
+            domain_yend = y_coords(yRows)              ! Last y value in first x-slice
+            ! Calculate simulation grid steps
+            x_step = x_cc(1) - x_cc(0)
+            y_step = y_cc(1) - y_cc(0)
+            ! Calculate offsets
+            delta_x = x_cc(index_x) - domain_xstart + x_step/2.0_wp
+            delta_y = y_cc(index_y) - domain_ystart + y_step/2.0_wp
+            ! Global offsets in case of multiple ranks
+            global_offset_x = nint(abs(delta_x)/x_step)
+            global_offset_y = nint(abs(delta_y)/y_step)
+
+            files_loaded = .true.
+        end if
+
+        idx = i + 1 + global_offset_x - index_x  ! x-direction index in data grid
+        idy = j + 1 + global_offset_y - index_y  ! y-direction index in data grid
+
+        do f = 1, sys_size - 1
+            if (f >= momxe) then
+                jump = 1
+            else
+                jump = 0
+            end if
+            q_prim_vf(f + jump)%sf(i, j, k) = stored_values(idx, idy, f)
+        end do
+        ! Set z velocity to zero
+        q_prim_vf(momxe)%sf(i, j, k) = 0.0_wp
+
         ! Put your variable assignments here
     case default
         call s_int_to_str(patch_id, iStr)

src/pre_process/include/ExtrusionHardcodedIC.fpp

@@ -0,0 +1,45 @@
+!> @brief Deallocate and reset variables used for IC extrusion.
+!>
+!> @details
+!>   Deallocates:
+!>     - stored_values(:,:,:).
+!>     - x_coords(:).
+!>     - y_coords(:).
+!>
+!>   Resets:
+!>     - files_loaded => .false.
+!>     - index_x, index_y, global_offset_x => 0.
+!>
+!>   (Note: The corresponding loader macro reads
+!>    `prim.<variable>.00.<timestep>.dat` files—pattern controlled by `zeros_default`—
+!>    from `examples/Case_File/D`.)
+
+#:def HardcodedDimensionsExtrusion()
+    integer :: xRows, yRows, nRows, iix, iiy
+    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/YourDirectory" ! 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 = "00000"  ! Default zeros (can be changed)
+#: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