add elliptic smoothing

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 |

src/common/m_variables_conversion.fpp

@@ -746,14 +746,14 @@ contains
 
         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 @@ contains
     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 @@ module m_global_parameters
     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 @@ module m_global_parameters
 
     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 @@ module m_global_parameters
     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,12 @@ module m_global_parameters
     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 +341,8 @@ contains
         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 +353,8 @@ contains
         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 +514,8 @@ contains
 
         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 +798,32 @@ contains
         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 @@ contains
         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 @@ contains
         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 @@ contains
         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)