add perturbed IC stuff

Author: wilfonba

src/post_process/m_mpi_proxy.fpp

@@ -110,6 +110,7 @@ contains
             & '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
 
         call MPI_BCAST(flux_wrt(1), 3, MPI_LOGICAL, 0, MPI_COMM_WORLD, ierr)

src/pre_process/include/3dHardcodedIC.fpp

@@ -7,6 +7,45 @@
     ! Case 302 - IGR Jet
     real(wp) :: r, ux_th, ux_am, p_th, p_am, rho_th, rho_am, y_th, z_th, r_th, eps_smooth
 
+    real(wp), allocatable, dimension(:, :) :: ih
+    integer :: i, j, pos, start, end
+    logical :: file_exist
+    character(len=10000) :: line
+    character(len=25) :: value
+
+    if (patch_icpp(patch_id)%hcid == 303) then
+        allocate(ih(0:m_glb, 0:p_glb))
+
+        if (interface_file == '.') then
+            call s_mpi_abort("Error: interface_file must be specified for hcid=303")
+        else
+            inquire (file=trim(interface_file), exist=file_exist)
+            if (file_exist) then
+                open(unit=10, file=trim(interface_file), status="old", action="read")
+                do i = 0, m_glb
+                    read(10, '(A)') line  ! Read a full line as a string
+                    start = 1
+
+                    do j = 0, p_glb
+                        end = index(line(start:), ',')  ! Find the next comma
+                        if (end == 0) then
+                            value = trim(adjustl(line(start:)))  ! Last value in the line
+                        else
+                            value = trim(adjustl(line(start:start+end-2)))  ! Extract substring
+                            start = start + end  ! Move to next value
+                        end if
+                        read(value, *) ih(i, j)  ! Convert string to numeric value
+                        if (.not. f_is_default(normMag)) ih(i,j )= ih(i,j) * normMag
+                        if (.not. f_is_default(normFac)) ih(i,j) = ih(i,j) + normFac
+                    end do
+                end do
+                close(10)
+            else
+                call s_mpi_abort("Error: interface_file specified for hcid=303 does not exist")
+            end if
+        end if
+    end if
+
     eps = 1e-9_wp
 
 #:enddef
@@ -86,6 +125,23 @@
 
         q_prim_vf(E_idx)%sf(i, j, k) = p_th*f_cut_on(r - r_th, eps_smooth)*f_cut_on(x_cc(i), eps_smooth) + p_am
 
+    case (303) ! 3D Interface from file
+
+        alph = 0.5_wp * (1 + (1._wp - 2._wp * eps) * &
+                    tanh((ih(start_idx(1) + i,start_idx(3) + k)  - y_cc(j))*0.1_wp))
+
+        q_prim_vf(advxb)%sf(i,j,k) = alph
+        q_prim_vf(advxe)%sf(i,j,k) = 1._wp - alph
+
+        q_prim_vf(contxb)%sf(i,j,k) = q_prim_vf(advxb)%sf(i,j,k) * 1._wp
+        q_prim_vf(contxe)%sf(i,j,k) = q_prim_vf(advxe)%sf(i,j,k) * (1._wp / 950._wp)
+
+        q_prim_vf(E_idx)%sf(i,j,k) = p0 + &
+            (q_prim_vf(contxb)%sf(i,j,k) + q_prim_vf(contxe)%sf(i,j,k)) * g0 * &
+            (ih(start_idx(1) + i, start_idx(3) + k) - y_cc(j))
+
+        if (surface_tension) q_prim_vf(c_idx)%sf(i,j,k) = alph
+
     case (370)
         ! This hardcoded case extrudes a 2D profile to initialize a 3D simulation domain
         @: HardcodedReadValues()

src/pre_process/m_mpi_proxy.fpp

@@ -72,7 +72,7 @@ contains
             & 'perturb_flow_mag', 'pref', 'rhoref', 'poly_sigma', 'R0ref',     &
             & 'Web', 'Ca', 'Re_inv', 'sigR', 'sigV', 'rhoRV', 'palpha_eps',    &
             & 'ptgalpha_eps', 'sigma', 'pi_fac', 'mixlayer_vel_coef', 'Bx0',   &
-            & 'mixlayer_perturb_k0']
+            & 'mixlayer_perturb_k0', 'normFac', 'normMag', 'p0', 'g0']
             call MPI_BCAST(${VAR}$, 1, mpi_p, 0, MPI_COMM_WORLD, ierr)
         #:endfor