Heat equation smoothing of IC

.fortitude.toml

@@ -1,5 +1,5 @@
 [check]
 # Ignored Rules and justification:
-ignore = ["E001","S001","S101","M011","F001","S041","T001","S101","S102","T002","T011","T042","M001"]
+ignore = ["E001","S001","S101","C121","S091","MOD011","C001","S101","S102","C002","PORT011","C072","C003","C131","C141","C092"]
 file-extensions = ["f90","fpp","fypp"]
 output-format = "pylint"

docs/documentation/case.md

@@ -269,6 +269,11 @@ These physical parameters must be consistent with fluid material's parameters de
 
 - `model_spc` and `model_threshold` are ray-tracing parameters. `model_spc` defines the number of rays per cell to render the model. `model_threshold` defines the ray-tracing threshold at which the cell is marked as the model.
 
+#### Elliptic Smoothing
+
+Initial conditions in which not all patches support the `patch_icpp(j)%smoothen` parameter can still be smoothed by applying iterations of the heat equation to the initial condition.
+This is enabled by adding `'elliptic_smoothing': "T",` and `'elliptic_smoothing_iters': N,` to the case dictionary, where `N` is the number of smoothing iterations to apply.
+
 ### 4. Immersed Boundary Patches
 
 | Parameter            | Type    | Description |

examples/1D_inert_shocktube/README.md

@@ -1,7 +1,7 @@
 # 1D Multi-Component Inert Shock Tube
 
 Reference: 
-> P. J. Martínez Ferrer, R. Buttay, G. Lehnasch, and A. Mura, “A detailed verification procedure for compressible reactive multicomponent Navier–Stokes solvers”, Comput. & Fluids, vol. 89, pp. 88–110, Jan. 2014. Accessed: Oct. 13, 2024. [Online]. Available: https://doi.org/10.1016/j.compfluid.2013.10.014
+> P. J. Martínez Ferrer, R. Buttay, G. Lehnasch, and A. Mura, “A detailed verification procedure for compressible reactive multicomponent Navier–Stokes solvers”, Computers & Fluids, vol. 89, pp. 88–110, Jan. 2014. Accessed: Oct. 13, 2024. [Online]. Available: https://doi.org/10.1016/j.compfluid.2013.10.014
 
 ## Initial Condition
 

examples/1D_reactive_shocktube/README.md

@@ -1,7 +1,7 @@
 # 1D Multi-Component Reactive Shock Tube
 
 References:
-> P. J. Martínez Ferrer, R. Buttay, G. Lehnasch, and A. Mura, “A detailed verification procedure for compressible reactive multicomponent Navier–Stokes solvers”, Comput. & Fluids, vol. 89, pp. 88–110, Jan. 2014. Accessed: Oct. 13, 2024. [Online]. Available: https://doi.org/10.1016/j.compfluid.2013.10.014
+> P. J. Martínez Ferrer, R. Buttay, G. Lehnasch, and A. Mura, “A detailed verification procedure for compressible reactive multicomponent Navier–Stokes solvers”, Computers & Fluids, vol. 89, pp. 88–110, Jan. 2014. Accessed: Oct. 13, 2024. [Online]. Available: https://doi.org/10.1016/j.compfluid.2013.10.014
 
 > H. Chen, C. Si, Y. Wu, H. Hu, and Y. Zhu, “Numerical investigation of the effect of equivalence ratio on the propagation characteristics and performance of rotating detonation engine”, Int. J. Hydrogen Energy, Mar. 2023. Accessed: Oct. 13, 2024. [Online]. Available: https://doi.org/10.1016/j.ijhydene.2023.03.190
 

misc/acc_devices.f90

@@ -23,4 +23,4 @@ program main
         print '(" - "I3" : "A" | "F0.2" GB | "A"")', i, devvendor, devmem, devname
     end do
 
-end program
+end program main

src/common/m_nvtx.f90

@@ -32,15 +32,15 @@ module m_nvtx
         subroutine nvtxRangePushA(name) bind(C, name='nvtxRangePushA')
             use iso_c_binding
 
-            character(kind=c_char, len=*) :: name
+            character(kind=c_char, len=*), intent(IN) :: name
         end subroutine nvtxRangePushA
 
         ! push range with custom label and custom color
         subroutine nvtxRangePushEx(event) bind(C, name='nvtxRangePushEx')
             use iso_c_binding
 
             import :: nvtxEventAttributes
-            type(nvtxEventAttributes) :: event
+            type(nvtxEventAttributes), intent(IN) :: event
         end subroutine nvtxRangePushEx
     end interface nvtxRangePush
 
@@ -54,8 +54,8 @@ end subroutine nvtxRangePop
 contains
 
     subroutine nvtxStartRange(name, id)
-        character(kind=c_char, len=*) :: name
-        integer, optional :: id
+        character(kind=c_char, len=*), intent(IN) :: name
+        integer, intent(IN), optional :: id
         type(nvtxEventAttributes) :: event
 
 #if defined(MFC_OpenACC) && defined(__PGI)

src/common/m_variables_conversion.fpp

@@ -746,14 +746,14 @@
 
         type(scalar_field), dimension(sys_size), intent(in) :: qK_cons_vf
 
-        real(wp), dimension(idwbuff(1)%beg:, idwbuff(2)%beg:, idwbuff(3)%beg:, 1:, 1:), intent(inout) :: mv
+        real(wp), dimension(idwint(1)%beg:, idwint(2)%beg:, idwint(3)%beg:, 1:, 1:), intent(inout) :: mv
 
         integer :: i, j, k, l
         real(wp) :: mu, sig, nbub_sc
 
-        do l = idwbuff(3)%beg, idwbuff(3)%end
-            do k = idwbuff(2)%beg, idwbuff(2)%end
-                do j = idwbuff(1)%beg, idwbuff(1)%end
+        do l = idwint(3)%beg, idwint(3)%end
+            do k = idwint(2)%beg, idwint(2)%end
+                do j = idwint(1)%beg, idwint(1)%end
 
                     nbub_sc = qK_cons_vf(bubxb)%sf(j, k, l)
 
@@ -778,15 +778,15 @@
     subroutine s_initialize_pb(qK_cons_vf, mv, pb)
         type(scalar_field), dimension(sys_size), intent(in) :: qK_cons_vf
 
-        real(wp), dimension(idwbuff(1)%beg:, idwbuff(2)%beg:, idwbuff(3)%beg:, 1:, 1:), intent(in) :: mv
-        real(wp), dimension(idwbuff(1)%beg:, idwbuff(2)%beg:, idwbuff(3)%beg:, 1:, 1:), intent(inout) :: pb
+        real(wp), dimension(idwint(1)%beg:, idwint(2)%beg:, idwint(3)%beg:, 1:, 1:), intent(in) :: mv
+        real(wp), dimension(idwint(1)%beg:, idwint(2)%beg:, idwint(3)%beg:, 1:, 1:), intent(inout) :: pb
 
         integer :: i, j, k, l
         real(wp) :: mu, sig, nbub_sc
 
-        do l = idwbuff(3)%beg, idwbuff(3)%end
-            do k = idwbuff(2)%beg, idwbuff(2)%end
-                do j = idwbuff(1)%beg, idwbuff(1)%end
+        do l = idwint(3)%beg, idwint(3)%end
+            do k = idwint(2)%beg, idwint(2)%end
+                do j = idwint(1)%beg, idwint(1)%end
 
                     nbub_sc = qK_cons_vf(bubxb)%sf(j, k, l)
 

src/pre_process/m_checker.fpp

@@ -192,6 +192,9 @@ contains
         @:PROHIBIT(mixlayer_perturb .and. num_fluids > 1, "mixlayer_perturb requires num_fluids = 1")
         @:PROHIBIT(mixlayer_perturb .and. any((/bc_y%beg, bc_y%end/) /= -6), &
             "mixlayer_perturb requires both bc_y%beg and bc_y%end to be 6")
+        @:PROHIBIT(elliptic_smoothing .and. elliptic_smoothing_iters < 1, &
+            "elliptic_smoothing_iters must be positive")
+
     end subroutine s_check_inputs_misc
 
 end module m_checker

src/pre_process/m_global_parameters.fpp

@@ -82,6 +82,7 @@
     integer :: num_fluids            !< Number of different fluids present in the flow
     logical :: mpp_lim               !< Alpha limiter
     integer :: sys_size              !< Number of unknowns in the system of equations
+    integer :: weno_polyn     !< Degree of the WENO polynomials (polyn)
     integer :: weno_order            !< Order of accuracy for the WENO reconstruction
     logical :: hypoelasticity        !< activate hypoelasticity
     logical :: hyperelasticity       !< activate hyperelasticity
@@ -130,6 +131,9 @@
 
     real(wp) :: pi_fac !< Factor for artificial pi_inf
 
+    logical :: viscous
+    logical :: bubbles_lagrange
+
     ! Perturb density of surrounding air so as to break symmetry of grid
     logical :: perturb_flow
     integer :: perturb_flow_fluid   !< Fluid to be perturbed with perturb_flow flag
@@ -138,6 +142,9 @@
     integer :: perturb_sph_fluid    !< Fluid to be perturbed with perturb_sph flag
     real(wp), dimension(num_fluids_max) :: fluid_rho
 
+    logical :: elliptic_smoothing
+    integer :: elliptic_smoothing_iters
+
     integer, allocatable, dimension(:) :: proc_coords !<
     !! Processor coordinates in MPI_CART_COMM
 
@@ -248,6 +255,11 @@
     type(pres_field) :: pb
     type(pres_field) :: mv
 
+    integer :: buff_size !<
+    !! The number of cells that are necessary to be able to store enough boundary
+    !! conditions data to march the solution in the physical computational domain
+    !! to the next time-step.
+
 contains
 
     !>  Assigns default values to user inputs prior to reading
@@ -328,6 +340,8 @@
         parallel_io = .false.
         file_per_process = .false.
         precision = 2
+        viscous = .false.
+        bubbles_lagrange = .false.
         mixlayer_vel_profile = .false.
         mixlayer_vel_coef = 1._wp
         mixlayer_domain = 1._wp
@@ -338,6 +352,8 @@
         perturb_sph = .false.
         perturb_sph_fluid = dflt_int
         fluid_rho = dflt_real
+        elliptic_smoothing_iters = dflt_int
+        elliptic_smoothing = .false.
 
         ! Initial condition parameters
         num_patches = dflt_int
@@ -497,6 +513,8 @@
 
         integer :: i, j, fac
 
+        weno_polyn = (weno_order - 1)/2
+
         ! Determining the layout of the state vectors and overall size of
         ! the system of equations, given the dimensionality and choice of
         ! the equations of motion
@@ -779,12 +797,32 @@
         chemxb = species_idx%beg
         chemxe = species_idx%end
 
+        ! Determining the number of cells that are needed in order to store
+        ! sufficient boundary conditions data as to iterate the solution in
+        ! the physical computational domain from one time-step iteration to
+        ! the next one
+        if (viscous) then
+            buff_size = 2*weno_polyn + 2
+        else
+            buff_size = weno_polyn + 2
+        end if
+
+        ! Correction for smearing function in the lagrangian subgrid bubble model
+        if (bubbles_lagrange) then
+            buff_size = max(buff_size, 6)
+        end if
+
         ! Configuring Coordinate Direction Indexes
         idwint(1)%beg = 0; idwint(2)%beg = 0; idwint(3)%beg = 0
         idwint(1)%end = m; idwint(2)%end = n; idwint(3)%end = p
 
-        ! There is no buffer region in pre_process.
-        idwbuff(1) = idwint(1); idwbuff(2) = idwint(2); idwbuff(3) = idwint(3)
+        idwbuff(1)%beg = -buff_size
+        if (num_dims > 1) then; idwbuff(2)%beg = -buff_size; else; idwbuff(2)%beg = 0; end if
+        if (num_dims > 2) then; idwbuff(3)%beg = -buff_size; else; idwbuff(3)%beg = 0; end if
+
+        idwbuff(1)%end = idwint(1)%end - idwbuff(1)%beg
+        idwbuff(2)%end = idwint(2)%end - idwbuff(2)%beg
+        idwbuff(3)%end = idwint(3)%end - idwbuff(3)%beg
 
 #ifdef MFC_MPI
 

src/pre_process/m_initial_condition.fpp

@@ -75,7 +75,9 @@
         allocate (q_cons_vf(1:sys_size))
 
         do i = 1, sys_size
-            allocate (q_prim_vf(i)%sf(0:m, 0:n, 0:p))
+            allocate (q_prim_vf(i)%sf(idwbuff(1)%beg:idwbuff(1)%end, &
+                                      idwbuff(2)%beg:idwbuff(2)%end, &
+                                      idwbuff(3)%beg:idwbuff(3)%end))
             allocate (q_cons_vf(i)%sf(0:m, 0:n, 0:p))
         end do
 
@@ -133,15 +135,15 @@
         character(len=10) :: iStr
 
         ! First, compute the temperature field from the conservative variables.
-        if (chemistry) call s_compute_q_T_sf(q_T_sf, q_cons_vf, idwbuff)
+        if (chemistry) call s_compute_q_T_sf(q_T_sf, q_cons_vf, idwint)
 
         ! Converting the conservative variables to the primitive ones given
         ! preexisting initial condition data files were read in on start-up
         if (old_ic) then
             call s_convert_conservative_to_primitive_variables(q_cons_vf, &
                                                                q_T_sf, &
                                                                q_prim_vf, &
-                                                               idwbuff)
+                                                               idwint)
         end if
 
         !  3D Patch Geometries
@@ -337,6 +339,7 @@
         if (perturb_flow) call s_perturb_surrounding_flow(q_prim_vf)
         if (perturb_sph) call s_perturb_sphere(q_prim_vf)
         if (mixlayer_perturb) call s_superposition_instability_wave(q_prim_vf)
+        if (elliptic_smoothing) call s_elliptic_smoothing(q_prim_vf)
 
         ! Converting the primitive variables to the conservative ones
         call s_convert_primitive_to_conservative_variables(q_prim_vf, q_cons_vf)